我在另一個平臺（這里）上問了一個問題——如果能得到您的意見，讓我的 Python 代碼在很短的時間內運行，那就太好了。目前，一個包含數百萬條目的檔案需要 3 個多小時。
from Bio import SeqIO
import sys
def QIAseq_UMI_correction():
script=sys.argv[0]
file_name=sys.argv[1]
dicts1 = {}
dicts2 = {}
lst = []
with open(file_name, "r") as Fastq:
for record in SeqIO.parse(Fastq,'fastq'):
#print(record.id)
#print(record.seq)
#looking for the 3 prime adapter
if "AACTGTAGGCACCATCAAT" in record.seq:
adapter_pos = record.seq.find('AACTGTAGGCACCATCAAT')
#Only record is used to be able to save the all atributes like phred score in the fastq file
miRNAseq = record[:adapter_pos]
adapter_seq=record[adapter_pos:adapter_pos 19]
umi_seq = record[adapter_pos 19:adapter_pos 19 12]
i = record.id
x = miRNAseq.seq umi_seq.seq
#print((miRNAseq umi_seq).format("fastq"))
dicts1[i]=x
#write ids and seq in a dictionary and keep one if there are multiple seqs with miRNA-seq UMI
for k,v in dicts1.items():
if v not in dicts2.values():
dicts2[k] = v
#making a list
for keys in dicts2:
lst.append(keys)
#print(dicts1)
#print(dicts2)
#print(lst)
with open(file_name, "r") as Fastq:
for record in SeqIO.parse(Fastq,'fastq'):
#based on the saved ids in the list print the entries (miRNA 3' adapter UMI)
if record.id in lst:
adapter_pos = record.seq.find('AACTGTAGGCACCATCAAT')
miRNAseq = record[:adapter_pos]
adapter_seq=record[adapter_pos:adapter_pos 19]
umi_seq = record[adapter_pos 19:adapter_pos 19 12]
#print(record.seq)
#print(miRNAseq.seq)
#print(adapter_seq.seq)
#print(umi_seq.seq)
#print("@" record.id)
if len(miRNAseq.seq adapter_seq.seq umi_seq.seq) <= 50:
print((miRNAseq adapter_seq umi_seq).format("fastq"),end='')
if len(miRNAseq.seq adapter_seq.seq umi_seq.seq) > 50:
cut = len(miRNAseq.seq adapter_seq.seq umi_seq.seq) - 50
print((miRNAseq adapter_seq umi_seq)[:-cut].format("fastq"), end='')
if __name__ == '__main__':
QIAseq_UMI_correction()
uj5u.com熱心網友回復：
為什么不對檔案進行一次讀取、決議和回圈呢？我已將第二個回圈的代碼移至第一個回圈，我是否遺漏了什么？為什么要回圈兩次？
from Bio import SeqIO
import sys
def QIAseq_UMI_correction():
script=sys.argv[0]
file_name=sys.argv[1]
dicts1 = {}
dicts2 = {}
lst = []
sentinel = 100
with open(file_name, "r") as Fastq:
for record in SeqIO.parse(Fastq,'fastq'):
# only for testing
if sentinel < 0:
break
sentinel -= 1
#print(record.id)
#print(record.seq)
#looking for the 3 prime adapter
if "AACTGTAGGCACCATCAAT" in record.seq:
adapter_pos = record.seq.find('AACTGTAGGCACCATCAAT')
#Only record is used to be able to save the all atributes like phred score in the fastq file
miRNAseq = record[:adapter_pos]
adapter_seq=record[adapter_pos:adapter_pos 19]
umi_seq = record[adapter_pos 19:adapter_pos 19 12]
i = record.id
x = miRNAseq.seq umi_seq.seq
#print((miRNAseq umi_seq).format("fastq"))
if x not in dicts2:
if len(miRNAseq.seq adapter_seq.seq umi_seq.seq) <= 50:
print((miRNAseq adapter_seq umi_seq).format("fastq"),end='')
if len(miRNAseq.seq adapter_seq.seq umi_seq.seq) > 50:
cut = len(miRNAseq.seq adapter_seq.seq umi_seq.seq) - 50
print((miRNAseq adapter_seq umi_seq)[:-cut].format("fastq"), end='')
dicts1[i]=x
dicts2[x]=i
if __name__ == '__main__':
QIAseq_UMI_correction()
其他建議：
如評論中所述，您可以對主要步驟進行計時，以查看可以縮短時間的地方。以timeit為例。我的建議是計時SeqIO.parse，if "AACTGTAGGCACCATCAAT" in record.seq:
我懷疑大部分時間都花在使用 SeqIO.parse 進行決議上，所以理想情況下你應該使用一次。
最后的建議是使用較小的記錄集，直到您的代碼準備好您需要它做的事情。我添加了一個哨兵變數作為示例，以便在探索 100 條匹配記錄時跳出回圈。
uj5u.com熱心網友回復：
我第一次看到的東西是：
首先檢查的地方
if "AACTGTAGGCACCATCAAT" in record.seq:
adapter_pos = record.seq.find('AACTGTAGGCACCATCAAT')
您可以使用以下內容來避免兩次搜索序列。
adapter_pos = record.seq.find('AACTGTAGGCACCATCAAT')
if adapter_pos != -1: # check that sequence is found
這些行：
i = record.id
x = miRNAseq.seq umi_seq.seq
#print((miRNAseq umi_seq).format("fastq"))
dicts1[i]=x
可以改為：
dicts1[record.id]=miRNAseq.seq umi_seq.seq
接下來的幾行：
for k,v in dicts1.items():
if v not in dicts2.values():
dicts2[k] = v
可以改為
dict2 = {**dict1,**dict2}
但是，我不確定這種變化實際上可能會以性能為代價。
接下來是那個
for keys in dicts2:
lst.append(keys)
可以洗掉和
lst = list(dicts2.keys())
可以在外部和第一次通過之后添加（您有注釋掉的列印陳述句。）
最后，正如@Roeften 建議的那樣，您可以將第二塊代碼放在第一塊代碼中，以避免兩次遍歷整個檔案。
至少其中一些建議會有所幫助，但實際上 python 并不是一種快速的語言，如果您想定期進行此類分析，您可能會考慮使用更快的方法。
uj5u.com熱心網友回復：
我會在這里使用磁區方法來回答您自己的問題。
from Bio import SeqIO
import sys
def QIAseq_UMI_correction():
script=sys.argv[0]
file_name=sys.argv[1]
dicts = {}
lst = []
with open(file_name, "r") as Fastq:
for record in SeqIO.parse(Fastq,'fastq'):
miRNAseq, adapter_seq, umi_seq = str(record.seq).partition("AACTGTAGGCACCATCAAT")
if adapter_seq:
x = miRNAseq.seq umi_seq.seq
y = miRNAseq.seq adapter_seq.seq umi_seq.seq
#print((miRNAseq umi_seq).format("fastq"))
dicts[miRNAseq umi_seq] = record.id
if len(y) <= 50:
print((miRNAseq adapter_seq umi_seq).format("fastq"),end='')
else:
cut = len(y) - 50
print((miRNAseq adapter_seq umi_seq)[:-cut].format("fastq"), end='')
if __name__ == '__main__':
QIAseq_UMI_correction()
uj5u.com熱心網友回復：
我通過洗掉第二部分來解決它 - 然后檢查 x 中是否不存在dicts.keys()（搜索dicts.values()速度要慢得多并且是真正的瓶頸）然后寫入y輸出。我現在只使用一個字典，代碼已經在生成輸出。這是更新的版本。
from Bio import SeqIO
import sys
def QIAseq_UMI_correction():
script=sys.argv[0]
file_name=sys.argv[1]
dicts = {}
lst = []
with open(file_name, "r") as Fastq:
for record in SeqIO.parse(Fastq,'fastq'):
if "AACTGTAGGCACCATCAAT" in record.seq:
adapter_pos = record.seq.find('AACTGTAGGCACCATCAAT')
miRNAseq = record[:adapter_pos]
adapter_seq=record[adapter_pos:adapter_pos 19]
umi_seq = record[adapter_pos 19:adapter_pos 19 12]
i = record.id
x = miRNAseq.seq umi_seq.seq
y = miRNAseq.seq adapter_seq.seq umi_seq.seq
#print((miRNAseq umi_seq).format("fastq"))
if x not in dicts.keys():
dicts[x]=i
if len(y) <= 50:
print((miRNAseq adapter_seq umi_seq).format("fastq"),end='')
if len(y) > 50:
cut = len(y) - 50
print((miRNAseq adapter_seq umi_seq)[:-cut].format("fastq"), end='')
if __name__ == '__main__':
QIAseq_UMI_correction()
轉載請註明出處，本文鏈接：https://www.uj5u.com/shujuku/432700.html
下一篇：診斷postgres中缺少的索引