SILVA、GREENGENES、RDP三大數據庫的序列探索統計
- 2020 年 3 月 3 日
- 筆記
最近對16s的三大數據庫的序列的具體序列情況挺好奇的，決定統計一下各個序列的長度分佈情況，以及這些序列具體分佈在哪幾個V區，有助於我解決後面16So數據的問題。還是用上我三腳貓的功夫，開始今天 的探索，沒有人探索的事情，還是挺開心的。
1.統計序列長度分佈情況
01 #獲得長度列表文件 02 length_list = [] 03 with open('current_Bacteria_unaligned.fa') as f: 04 flag = 0 05 length = 0 06 for line in f: 07 if line.startswith('>'): 08 flag = 1 09 if length != 0: 10 #print(length) 11 length_list.append(length) 12 #break 13 length = 0 14 else: 15 length = 0 16 elif flag == 1: 17 length +=len(line) 18 19 fout = open('length_table.txt', 'w') 20 for a in length_list: 21 fout.write(str(a) + 'n') 22 fout.close() 01 #統計長度區間分佈 02 length_150 = 0 03 length_300 = 0 04 length_600 = 0 05 length_1300 = 0 06 length_1300_ = 0 07 with open('length_table.txt') as f: 08 for line in f: 09 if 0 < int(line.strip()) <= 150: 10 length_150 += 1 11 elif 150 < int(line.strip()) <= 300: 12 length_300 += 1 13 elif 300 < int(line.strip()) <= 600: 14 length_600 += 1 15 elif 600 < int(line.strip()) <= 1300: 16 length_1300 += 1 17 elif int(line.strip()) > 1300: 18 length_1300_ += 1 19 print('150以內：', length_150,'n', 20 '150-300:', length_300, 'n', 21 '300-600:', length_600, 'n', 22 '600-1300:', length_1300,'n', 23 '1300-:', length_1300_) 1 #最後結果 2 150以內： 0 3 150-300: 0 4 300-600: 617,554 5 600-1300: 1,063,377 6 1300-: 1,515,109
我也可視化一下：
更正一下，這裡用的是RDP數據庫，之前由於Silva數據庫用的是不兼容的14級分類系統而沒採用。
從圖中可以看出，大部分序列還是集中在600以上。
接着是greengenes數據庫，這個數據庫雖然序列較少，但是長度大部分集中在1300+，質量較高，就是好久沒更新過了。 150以內： 0 150-300: 0 300-600: 0 600-1300: 895 1300-: 1262090
2.統計V區分佈情況
從一個公眾號得到的一張分佈圖是這樣的，
我想確定的是序列都包含在哪個或者哪兩個區。花了好大有力氣，才把代碼理好，效率應該還是比較低的，但是，對於我來說
001 import re 002 003 dic ={} 004 v1f = re.compile('GGATCCAGACTTTGATYMTGGCTCAG', re.I) 005 v3f = re.compile('CCTA[CT]GGG[AG][GTC]GCA[CG]CAG', re.I) 006 v4f = re.compile('GTG[CT]CAGC[AC]GCCGCGGTAA', re.I) 007 v4r = re.compile('ATTAGA[AT]ACCC[CTG][ATGC]GTAGTCC', re.I) 008 v6f = re.compile('AACGCGAAGAACCTTAC', re.I) 009 v8f = re.compile('CGTCATCC[AC]CACCTTCCTC', re.I) 010 vr = re.compile('AAGTCGTAACAAGGTA[AG]CCGTA', re.I) 011 #v4r = re.compile('GGACTAC[ATGC][ACG]GGGT[AT]TCTAAT', re.I) 012 vf = re.compile('AG[AG]GTT[CT]GAT[CT][AC]TGGCTCAG', re.I) 013 #vr = re.compile('GACGGGCGGTG[AT]GT[AG]CA', re.I) 014 #seq_list = [] 015 016 017 def decide_which_zone(f1, f3, f4, f4r, f6, f8, f, fr): 018 type = [f1, f3, f4, f4r, f6, f8, f, fr] 019 type_str = ['f1', 'f3', 'f4', 'f4r', 'f6', 'f8', 'f', 'fr'] 020 type_name = [] 021 for i in range(len(type)): 022 #print(type[i]) 023 if type[i] != None : 024 type_name.append(type_str[i]) 025 else: 026 continue 027 if i + 1 < len(type): 028 if type[i + 1] != None and type_str[i + 1] not in type_str: 029 type_name.append(type_str[i + 1]) 030 else: 031 continue 032 if i + 2 < len(type): 033 if type[i + 2] != None and type_str[i + 2] not in type_str: 034 type_name.append(type_str[i + 2]) 035 else: 036 continue 037 if i + 3 < len(type): 038 if type[i + 3] != None and type_str[i + 3] not in type_str: 039 type_name.append(type_str[i + 3]) 040 else: 041 continue 042 if i + 4 < len(type): 043 if type[i + 4] != None and type_str[i + 4] not in type_str: 044 type_name.append(type_str[i + 4]) 045 else: 046 continue 047 if i + 5 < len(type): 048 if type[i + 5] != None and type_str[i +5 ] not in type_str: 049 type_name.append(type_str[i + 5]) 050 else: 051 continue 052 if i + 6 < len(type): 053 if type[i + 6] != None and type_str[i + 6] not in type_str: 054 type_name.append(type_str[i + 6]) 055 else: 056 continue 057 if i + 7 < len(type): 058 if type[i + 7] != None and type_str[i + 7] not in type_str: 059 type_name.append(type_str[i + 7]) 060 else: 061 continue 062 return type_name 063 064 with open('current_Bacteria_unaligned.fa') as f: 065 li = ['f1', 'f3', 'f4', 'f4r', 'f6', 'f8', 'f', 'fr'] 066 for i in range(len(li)): 067 for j in range(len(li)): 068 if i <= j: 069 dic[li[i] + li[j]] = 0 070 071 072 #print(dic) 073 flag = 0 074 seq = '' 075 #i = 0 076 for line in f: 077 if line.startswith('>'):# and i <=2: 078 flag = 1 079 if seq != '': 080 #print(seq) 081 f1 = v1f.search(seq) 082 f4r = v4r.search(seq) 083 f3 = v3f.search(seq) 084 f4 = v4f.search(seq) 085 f6 = v6f.search(seq) 086 f8 = v8f.search(seq) 087 f = vf.search(seq) 088 fr = vr.search(seq) 089 #seq_list.append(length) 090 type_name = [] 091 type_name = decide_which_zone(f1, f3, f4, f4r, f6, f8, f, fr) 092 print(type_name) 093 if type_name != []: 094 type_key = type_name[0] + type_name[-1] 095 dic[type_key] += 1 096 #i += 1 097 #break 098 seq = '' 099 flag = 0 100 else: 101 seq = '' 102 elif flag == 1: 103 seq += line.strip() 104 print(dic) 105 #運行結果 106 {'f1f1': 0, 'f1f3': 0, 'f1f4': 1, 'f1f4r': 0, 'f1f6': 0, 'f1f8': 0, 'f1f': 0, 'f1fr': 0, 'f3f3': 120867, 'f3f4': 143686, 'f3f4r': 356027, 'f3f6': 400108, 'f3f8': 4, 'f3f': 171637, 'f3fr': 53716, 'f4f4': 17211, 'f4f4r': 71719, 'f4f6': 40223, 'f4f8': 4, 'f4f': 13185, 'f4fr': 3646, 'f4rf4r': 39537, 'f4rf6': 45160, 'f4rf8': 0, 'f4rf': 977, 'f4rfr': 4263, 'f6f6': 43983, 'f6f8': 0, 'f6f': 64, 'f6fr': 2300, 'f8f8': 0, 'f8f': 0, 'f8fr': 0, 'ff': 1670, 'ffr': 26, 'frfr': 3678}
上圖依然是RDP數據庫，前面標錯了。應該是某個地方出了問題，序列不全，少了估計有一半。不過趨勢還是比較明顯的，V3-V4是最多的，這也是前些年主流的研究的測序方式。
