2017-02-21 94 views
4

我試圖按照文件中序列的字母順序(而不是序列的ID)對fasta文件進行排序。 fasta文件包含200多個序列,我試圖在一個位主(使用python代碼)內找到重複數據(通過重複數據表示幾乎相同的蛋白質序列,但不是相同的ID)。 所以我想從fasta文件中創建一個字典,然後對字典的值進行排序。 我想使用的代碼如下:「NotImplementedError:SeqRecord」當使用排序在使用SeqIO解析的fasta文件時使用SeqIO

from Bio import SeqIO 


input_file = open("PP_Seq.fasta")  
my_dict = SeqIO.to_dict(SeqIO.parse(input_file, "fasta")) 
print sorted(my_dict.values()) 

我不斷收到此消息錯誤:

"Traceback (most recent call last): 
    File "sort.py", line 4, in <module> 
    print sorted(my_dict.values()) 
    File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/Bio/SeqRecord.py", line 730, in __lt__ 
    raise NotImplementedError(_NO_SEQRECORD_COMPARISON) 
NotImplementedError: SeqRecord comparison is deliberately not implemented. Explicitly compare the attributes of interest." 

我也試圖尋找如何鰭這個錯誤,但有多少ares't有關這方面的信息,以及我讀到的很少的信息,這些信息顯然是說存儲在詞典字典中的序列長度可能是一個問題?如果是這樣,如何在沒有SeqIO的情況下對fasta文件進行排序?

+0

你的字典是否像'{fasta_header:sequence}'? –

+1

這意味着'SeqRecords'不是可比較的,所以它們不能被排序。你想用什麼關鍵字來排序?像'sorted(my_dict.values(),key = operator.attrgetter('seq'))''可能會工作。 – mata

+0

@mata可以說我有這個文件作爲輸入: > seq0 ABCWYXO > SEQ1 IJKLMNOP > SEQ2 BCDEFGH > SEQ3 ABCDEFG 我所要的輸出是一個文件安排是這樣的: > SEQ3 ABCDEFG > seq0 ABCWYXO > SEQ2 BCDEFGH > SEQ1 IJKLMNOP 基本上由蛋白質序列的字母順序排序他們了...所以我想比較喜歡一個序列(〜應變g)在一個循環中對另一個字符進行排序,並以這種方式對它們進行排序。能夠根據該順序將它們放置在新文件中,並且每次都檢索自己的ID ... –

回答

2

正如mata說,你需要的關鍵功能傳遞給sorted

from Bio import SeqIO 
import operator 
input_file = open("example.fasta")  
my_dict = SeqIO.to_dict(SeqIO.parse(input_file, "fasta")) 
for r in sorted(my_dict.values(), key=operator.attrgetter('seq')): 
    print r.id, str(r.seq) 

回報:現在

seq3 ABCDEFG 
seq0 ABCWYXO 
seq2 BCDEFGH 
seq1 IJKLMNOP 

,你要完成什麼。如果您已按字母順序對200個序列進行排序,則仍需手動掃描列表才能找到相近的重複項。這很容易出錯,所以最好還是寫一些代碼。

在計算機科學中,edit distance是一種通過計算將一個字符串轉換爲另一個字符串所需的最小操作數來量化不同的兩個字符串(例如單詞)彼此之間的差異的方法。

該算法有幾種實現可用。我們將採取從this answer

def levenshteinDistance(s1, s2): 
    if len(s1) > len(s2): 
     s1, s2 = s2, s1 

    distances = range(len(s1) + 1) 
    for i2, c2 in enumerate(s2): 
     distances_ = [i2+1] 
     for i1, c1 in enumerate(s1): 
      if c1 == c2: 
       distances_.append(distances[i1]) 
      else: 
       distances_.append(1 + min((distances[i1], distances[i1 + 1], distances_[-1]))) 
     distances = distances_ 
    return distances[-1] 

現在我們需要決定兩個序列可能有多不同(多少個插入/去合/替代)的閾值。然後,我們兩兩在我們的FASTA文件中每兩個序列進行比較:

from Bio import SeqIO 
from itertools import combinations 
input_file = open("example.fasta")  

treshold = 4 
records = SeqIO.parse(input_file, "fasta") 
for record1, record2 in combinations(records, 2): 
    edit_distance = levenshteinDistance(str(record1.seq), str(record2.seq)) 
    if edit_distance <= treshold: 
     print "{} and {} differ in {} characters".format(record1.id, record2.id, edit_distance) 

這給:

seq0 and seq3 differ in 4 characters 
seq2 and seq3 differ in 2 characters 
1

這裏是爲了消除基於bioawkfastq-tools(也awk和FASTA文件複製另一種方式uniq通常存在於任何類似UNIX的命令行環境中):

bioawk -c fastx '{print "@"$name"\n"$seq"\n+\n"$seq}' test.fasta \ 
    | fastq-sort -s \ 
    | bioawk -c fastx '{print $name"\t"$seq}' \ 
    | uniq -f 1 \ 
    | awk '{print ">"$1"\n"$2}' 

bioawkawk的修改版本,其有助於操縱一些常見的生物信息學格式。

第一行將數據轉換爲fastq格式,因爲這是fastq-sort可以使用的格式。 -c fastx選項告訴bioawk將數據解析爲fasta或fastq。這定義了可在命令中使用的$name$seq字段。我們使用$seq字段作爲虛擬品質來獲得有效的fastq格式。

第二行告訴fastq-sort(來自fastq-tools)按順序對數據進行排序(選項-s)。

第三行使用bioawk提取的名稱和順序,並把它們作爲兩個製表符分隔的領域,在每個記錄一行的相關信息。

第四行使用uniq在比較連續行時(如果我正確理解-f選項的含義,我剛發現),忽略忽略第一個字段的重複項。如果我對uniq的理解是正確的,則融合記錄的名稱將是同序列中第一條記錄的名稱。

第五行使用awk重新格式化製表符分隔的字段FASTA成。

我測試了我的一些數據的這種做法,它似乎工作。