从fasta文件中找到唯一的fastq文件的第一行和第一行[关闭]答案

【问题标题】：find unique first top and bottom lines of fastq file from fasta file [closed]从fasta文件中找到唯一的fastq文件的第一行和第一行[关闭]
【发布时间】：2013-08-02 19:06:49
【问题描述】：

我有 2 个文件，一个是 fasta 文件，另一个是 fastq 文件。我想获取fasta，读取每一行并在fastq 文件中搜索每一行并打印顶行和底行。这就是我所拥有的

fasta 文件

读取1

啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAC

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAG

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGA

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGG

for seq in `cat sequences`;do grep -A 2 -B 1 $seq FAP.1.txt;done

@DH1DQQN1:269:C1UKCACXX:1:1107:20386:6577 1:N:0:TTAGGC

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAC

+

CCCFFFFFHGHHHJIJHFDDDB173@8815BDDB###############

@DH1DQQN1:269:C1UKCACXX:1:1114:5718:53821 1:N:0:TTAGGC

啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊

+ ;@?DBD

@DH1DQQN1:269:C1UKCACXX:1:1209:10703:35361 1:N:0:TTAGGC

啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊

+

@@@FFFFFHGHHHGIJHFDDDDDBDD69@6B-707537BDDDB75@@85

@DH1DQQN1:269:C1UKCACXX:1:1210:18926:75163 1:N:0:TTAGGC

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAG

@CCFFFFFHHHHHJJJHFDDD@77BDDDDB077007@B###########

从这里我们可以看到AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA出现了两次，但我只想打印一次。我该怎么做？

最终输出文件

@DH1DQQN1:269:C1UKCACXX:1:1107:20386:6577 1:N:0:TTAGGC

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAC

+

CCCFFFFFHGHHHJIJHFDDDB173@8815BDDB###############

@DH1DQQN1:269:C1UKCACXX:1:1114:5718:53821 1:N:0:TTAGGC

啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊

+

;@?DBD

@DH1DQQN1:269:C1UKCACXX:1:1210:18926:75163 1:N:0:TTAGGC

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAG

+

@CCFFFFFHHHHHJJJHFDDD@77BDDDDB077007@B

【问题讨论】：

生物学家你怎么了.. 总是发布格式最差的问题！您需要阅读 stackoverflow.com/help/formatting 以了解格式在 Stackoverflow 上的工作原理，以及在提问时阅读 stackoverflow.com/help/asking。
对不起！下次不会那样了
在过去的几个月里，您问了 14 个问题，所有问题的格式都很糟糕，必须彻底检修（似乎大部分是我自己做的）。如果您希望社区能够帮助您，那么您至少可以正确地格式化您的问题。

标签： python unix awk grep fastq

【解决方案1】：

除非您有充分的理由自己这样做，否则请使用Biopython。

快速：

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAC
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAG
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGA
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGG

fastq（基于您的但不完全相同，因为您的输出格式错误）：

@DH1DQQN1:269:C1UKCACXX:1:1107:20386:6577 1:N:0:TTAGGC
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAC
+
CCCFFFFFHGHHHJIJHFDDDB173@8815BDDB###############
@DH1DQQN1:269:C1UKCACXX:1:1114:5718:53821 1:N:0:TTAGGC
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
+
CCCFFFFFHGHHHJIJHFDDDB173@8815BDDB###############
@DH1DQQN1:269:C1UKCACXX:1:1209:10703:35361 1:N:0:TTAGGC
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
+
@@@FFFFFHGHHHGIJHFDDDDDBDD69@6B-707537BDDDB75@@85
@DH1DQQN1:269:C1UKCACXX:1:1210:18926:75163 1:N:0:TTAGGC
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAG
+
@CCFFFFFHHHHHJJJHFDDD@77BDDDDB077007@B###########

代码：

from Bio import SeqIO

with open("fasta") as fh:
    fasta = fh.read().splitlines()

seen = set()

for record in SeqIO.parse(open('fastq'), 'fastq'):
    seq = str(record.seq)
    if seq in fasta and seq not in seen:
        seen.add(seq)
        print record.format('fastq')

EDIT：上面是按fastq文件的顺序打印记录，而不是fasta文件。如果顺序不重要，则应使用该方法。否则，您可以将记录添加到字典中，其中键是它们在 FASTA 文件中的索引，最后将它们全部打印出来，对字典进行排序：

from Bio import SeqIO
import sys

with open("fasta") as fh:
    fasta = fh.read().splitlines()

seen = set()
records = {}

for record in SeqIO.parse(open('fastq'), 'fastq'):
    seq = str(record.seq)
    if seq in fasta and seq not in seen:
        seen.add(seq)
        records[fasta.index(seq)] = record

for record in sorted(records):
    sys.stdout.write(records[record].format('fastq'))

（这里我也使用sys.stdout.write而不是print，以避免额外的换行符。）

【讨论】：

我只想在fastq文件中打印一次fasta文件的第一行@paulo Almeida
这不是您在“最终输出文件”中的内容。您应该在问题中有一个所需输出的样本（我不明白您所说的“在 fastq 文件中打印一次 fasta 文件的第一行”是什么意思）。在任何情况下，一旦你拥有 Biopython 对象中的记录，你就可以只打印你想要的。
在最终的输出文件中，我提到第一行只出现了一次！这意味着我在 fastq 文件中重复了相同序列的行，我想用底行和顶行打印第一行 @paulo Alimeida
这就是这段代码的作用；它只打印该行一次（已打印的每个序列都被存储，因此不会打印两次）。
@abh，我的原始代码没有按照 fasta 文件对序列进行排序，现在可以了。

【解决方案2】：

所以听起来你想要一个只有唯一序列的 fastq？

这是一种非常低效的方法，但它应该可以工作。它将您的 fastq 文件存储为列表，因此希望它不会太大。它只会抛出重复的序列，而不是质量分数或任何东西。

fastqFile = list(open(fastq))
out = []
output = open('output.fastq', 'at')

for lineNum, line in enumerate(fastqFile):
    if lineNum < 4:
        out.append(line)
        output.write(line)
    else:
        if line not in out and lineNum % 4 != 3:
            output.write(fastqFile[lineNum - 1])
            output.write(line)
            output.write(fastqFile[lineNum + 1])
            output.write(fastqFile[lineNum + 2])
            out.append(fastqFile[lineNum - 1])
            out.append(line)
            out.append(fastqFile[lineNum + 1])
            out.append(fastqFile[lineNum + 2])

【讨论】：

【解决方案3】：

我想我知道你想说什么，所以这是我的代码。根据要求，它只会采用第一次出现的 fasta 序列。这可能不是最好的方法，但我是 python 新手。

# open the file into a list
fasta = open('fasta1.fa', 'r').read().splitlines()
fastq = open('fastq1.fq', 'r').read().splitlines()

# get rid of headers
# if headers important, please indicate in example
fastaseq = [s for s in fasta if not any('>' in t for t in s)]

# get rid of whitespace
fastaseq = filter(None, fastaseq)
fastq = filter(None, fastq)

# new list
newfastq = []

# go through each item in your fasta list
# if it matches, get the line above and below
# put in the new list
for fa in fastaseq:
    if fa in fastq:
        ind = fastq.index(fa)
        printblock = fastq[ind-1:ind+2]
    elif fa not in fastq:
        printblock = []
    if printblock:
        newfastq.append(printblock)

# print everything to file  
with open('fastq2.fq', 'w') as f:
    for block in newfastq:
        for item in block:
            f.writelines(item + '\n')

【讨论】：