生信中的Python——定位polyN
--這個轉(zhuǎn)自我的博客，我覺得很有意義，同時試一下B站的排版, https://atlasbioinfo.github.io/blog/--
馬上快畢業(yè)了，師妹跟我聊說她現(xiàn)在程序還是不太會寫，讓我在撤退之前給她補(bǔ)救一下，所以我準(zhǔn)備推出一系列生信代碼教程。主要是通過一些示例，寫一下我寫生信程序的思路和經(jīng)驗，并不包括具體的語法。第一個就寫一個有關(guān)polyN定位的例子，示例來自于已經(jīng)工作的孟師妹問我的一個問題，在此與大家分享。
生信案例描述
最直接的解決方案
適用于更多情況
寫在前面：
我并不想把python的語法或者什么再重復(fù)一遍，沒有任何意義。現(xiàn)在描述python基本語法的教學(xué)網(wǎng)站很多，搜索一大把。一般如果語法忘了，大家可以直接搜，搜索引擎也是IDE嘛。
1. 生信案例描述
故事是這樣的：
已經(jīng)工作的孟師妹突然給我發(fā)消息，說有一個task她的程序跑的特別慢，由于公司需要效率，所以問我能不能優(yōu)化一下程序。
她有一個基因組的fasta文件，里面包含poly N序列，現(xiàn)在需要把這些序列定位出來生成一個類似gff的表。
fasta示例文件如下：
>chr1
AACGNNNNACGTAC
那么，就應(yīng)該生成
#label? Chr? ? ?Beg End
polyN? ?chr1? ? 5? ?8? ?
簡單吧，那么，游戲開始~
注意：
我非常建議你先自己想一下怎么解決這個問題，寫好代碼再繼續(xù)看下一節(jié)。這樣收獲更大。
同理，每一節(jié)我還會提出問題，建議還是先完成相應(yīng)的代碼，再看下一節(jié)。
2. 最直接的解決方案
最直接的方案就是：
一行一行把文件讀出來
遇到“>”記錄一下fasta的title，序列開始計算
把序列一個個讀出來，看是不是“N”，如果是記錄開始位置，直到不是N時輸出。
這樣就出現(xiàn)了下面的代碼：
首先先建立一個叫test.fa的文件保存序列，內(nèi)容為：
>chr1
AACGNNNNACGTAC
之后寫個python
print("#label\tChr\tBeg\tEnd")
f=open("./test.fa",'r') #開文件
lines=f.readlines() #按照行讀取，保存在名為lines的list里面
#lines[0]是大于號的那行，把大于號后面的內(nèi)容輸出,strip()去換行符
print("polyN\t"+lines[0][1:].strip(),end="")
#之后需要處理序列那行，設(shè)置一個小標(biāo)簽jump，防止重復(fù)輸出
#0代表false，>0為true。
#按理說賦值數(shù)字就行，為了方便理解，我先用true和false代替
jump=True
for i in range(len(lines[1])):? #把每一個堿基拉過來遍歷
? ? if (lines[1][i] == "N"): #如果是N，記錄位置
? ? ? ? if (jump): #第一次遇到N，為true，記錄位置；之后為false，不再記錄
? ? ? ? ? ? print("\t"+str(i+1),end="") #記錄開始位置
? ? ? ? ? ? jump=False #以后在上面的jump判斷中不再通過，不再重復(fù)輸出起始位置
? ? else:
? ? ? ? if (not jump): # 如果這個堿基不是N，并且已經(jīng)開始記錄N位置了(jump為False)
? ? ? ? ? ? ? ? ? ? ? ? # 那么，not jump為true，記錄一下結(jié)束位置。
? ? ? ? ? ? print("\t"+str(i))
? ? ? ? ? ? break? ? ? ? ? ? ? ??
就像預(yù)期，輸出和上面示例輸出一模一樣：
#label? Chr? ? ?Beg? ? ?End
polyN? ?chr1? ? 5? ? ? ?8
評價一下上面的代碼。
上面內(nèi)容有一個小技巧就是jump，因為遇到第一個N需要記錄起始位置，直到第一個非N位點(diǎn)記錄終止位置。
不過，這個代碼已經(jīng)知道了需要讀取的文件有2行，第一行為名字，第二行為序列，且只有一段polyN，所以相當(dāng)于針對于這個數(shù)據(jù)“定制的”代碼。
這段代碼可以說"泛化"極弱，換個fasta數(shù)據(jù)立刻報錯。不過我寫這個教程的時候并沒有把這個最簡單的程序這個省去，因為這也是我寫程序的思路之一。
Atlas編程哲學(xué)
不要浪費(fèi)時間在沒有用的考慮上，一切都是為了最快的獲得正確的結(jié)果。
如果某個task以后應(yīng)用的場景極低，可以認(rèn)為是一次性代碼，用完就扔，那么就忘記“泛化”，最快速度達(dá)到目的才是最重要的。
?也就是說：如果你需要在這個test.fa后面再加點(diǎn)序列，直接打開編輯器加就行了，不需要寫一個通用的文檔追加程序，完善所有功能，做成一個產(chǎn)品，或者自己開發(fā)個編輯器啥的......
?3. 適用于更多情況
上一個程序的弊端在于，只能處理2行的，只有一段polyN的序列。真正的polyN序列往往很多段，fasta文件也有很多行，很多染色體，于是上面的程序就用不了了。
3.1 支持多染色體和多段polyN
所以首先需要增加的功能就是：
支持多個染色體在一個fasta中的情況
支持多段polyN
print("#label\tChr\tBeg\tEnd") #打印title
f=open("./test.fa",'r') #開文件
lines=f.readlines() #按照行讀取，保存在名為lines的list里面
chrname="" #設(shè)個變量記錄染色體名字
for line in lines:?
? ? line=line.strip() #去除一下末尾的換行
? ? if (line[:1] == '>'): #如果line的第一個是">"
? ? ? ? chrname=line[1:] #記錄一下染色體名字
? ? else:
? ? ? ? beg=-1 #首先設(shè)置一個叫beg的標(biāo)簽，記錄起始坐標(biāo)
? ? ? ? for i in range(len(line)):? #把每一個堿基拉過來遍歷
? ? ? ? ? ? if (line[i] == "N"): #如果是N，記錄位置
? ? ? ? ? ? ? ? if (beg == -1):
? ? ? ? ? ? ? ? ? ? beg=i+1
? ? ? ? ? ? else: #如果不是N，且beg已經(jīng)開始記錄了，那么輸出前一個位置
? ? ? ? ? ? ? ? if (beg>0):
? ? ? ? ? ? ? ? ? ? print("polyN\t"+chrname+"\t"+str(beg)+"\t"+str(i))
? ? ? ? ? ? ? ? ? ? beg=-1? #最后，beg初始化，賦值-1進(jìn)入下一輪循環(huán)
上段代碼輸出為：
#label? Chr? ? ?Beg? ? ?End
polyN? ?chr1? ? 5? ? ? ?8
woo，和期望結(jié)果的一模一樣。
同時，如果我們弄得稍微復(fù)雜點(diǎn)，有2段polyN，上面這段代碼也是OK的。比如，test.fa改成：
>chr1
AACGNNNNACGTNNNNACGTA
則輸出：
#label? Chr? ? ?Beg? ? ?End
polyN? ?chr1? ? 5? ? ? ?8
polyN? ?chr1? ? 13? ? ? 16
多染色體也OK，test.fa修改成：
>chr1
AACGNNNNACGTNNNNACGTA
>chr2
AACGNNNNACGTNNNNACGTA
則輸出：
#label? Chr? ? ?Beg? ? ?End
polyN? ?chr1? ? 5? ? ? ?8
polyN? ?chr1? ? 13? ? ? 16
polyN? ?chr2? ? 5? ? ? ?8
polyN? ?chr2? ? 13? ? ? 16
不過，如果fasta序列是這樣的：
>chr1
AACGNNNNACGTNNNNACGTA
AACGNNNNACGTNNNNACGTA
額，輸出是：
#label? Chr? ? ?Beg? ? ?End
polyN? ?chr1? ? 5? ? ? ?8
polyN? ?chr1? ? 13? ? ? 16
polyN? ?chr1? ? 5? ? ? ?8
polyN? ?chr1? ? 13? ? ? 16
也就是程序目前只能處理單行fasta，無法多行，所以需要進(jìn)一步修改。
3.2 支持多行
多行，首先想到的是需要一個變量記錄每行的長度，定位的位置加上這個長度不就行了。
pos=0 # for循環(huán)前定義一個pos變量
? ? else:
? ? ? ? ...#省略代碼
? ? ? ? beg = pos + i + 1
? ? ? ? ...#省略代碼
? ? pos = pos + len(line) #最后一行加這個
輸出為：
#label? Chr? ? ?Beg? ? ?End
polyN? ?chr1? ? 5? ? ? ?8
polyN? ?chr1? ? 13? ? ? 16
polyN? ?chr1? ? 26? ? ? 29
polyN? ?chr1? ? 34? ? ? 37
這下正常了，不過別急，如果換行呢？比如：
>chr1
ACTAGCATCGTACNNN
NNNACGTACGTACGTA
輸出：
#label? Chr? ? ?Beg? ? ?End
polyN? ?chr1? ? 17? ? ? 19
第一行直接掠過去了（由于沒有非N作為阻斷），直接記錄第二行，位置倒是對的。
所以，現(xiàn)在我們有2個策略：
1. 把多行fasta首先轉(zhuǎn)換成單行fasta，直接用前面的代碼就行。
2. 想出一個多行的策略解決這種情況
那么，應(yīng)該用哪個方案呢？
Atlas編程哲學(xué)
在邊界內(nèi)效率最大化。
軟件工程師，為啥叫工程師而不是軟件科學(xué)家？因為工程師需要考慮到實際情況，而不是所有研究全部處于理論狀態(tài)。
對于寫代碼，需要考慮的就是能用的內(nèi)存、硬盤空間和CPU，以及代碼運(yùn)行起來的時間復(fù)雜度和空間復(fù)雜度。
也就是說：
如果我所需要處理的情況，就是上面那個示例文件，我就直接肉眼數(shù)出來就行，也別寫代碼了。
如果我需要處理較大的染色體，而內(nèi)存可以塞下，那么我就采用方案1，簡單粗暴，多行變一行，直接解決。
如果我內(nèi)存不夠，我就需要采用方案2，修改策略，單行處理。
對于其他情況，通常還需要考慮自己的硬件優(yōu)勢
如果處理很多基因組，而我的優(yōu)勢的多CPU，就需要對文件進(jìn)行拆分進(jìn)行多進(jìn)程處理；或者寫多線程進(jìn)行處理。
如果優(yōu)勢是大硬盤，劣勢是小內(nèi)存，那么要么一點(diǎn)一點(diǎn)處理，要么把內(nèi)存的東西先寫到硬盤，一步一步處理。
總的來說，寫程序并不是黑寫，需要像工程師那樣考慮實際情況進(jìn)行效率最大化的優(yōu)化，可能在不同場景需要不同的處理策略。當(dāng)然，如果硬件很好，大內(nèi)存，大硬盤，主頻高核還多的CPU，就可以任性，想怎么寫就怎么寫。
回到主題，我傾向于寫多行的策略，適用范圍廣還挑戰(zhàn)一下自己；不過如果平時著急出結(jié)果的代碼，我會用方案1盡快出結(jié)果。
其實，改起來很簡單，把beg往外圈放就好了：
print("#label\tChr\tBeg\tEnd")?
f=open("./test.fa",'r')?
lines=f.readlines()?
chrname=""?
pos=0
beg=-1 #把beg移動到這個位置
for line in lines:?
? ? line=line.strip()?
? ? if (line[:1] == '>'):
? ? ? ? chrname=line[1:]?
? ? ? ? beg = -1 # 這里還需要再初始化一下，重置一下beg和pos
? ? ? ? pos=0
? ? else:
? ? ? ? for i in range(len(line)):?
? ? ? ? ? ? if (line[i] == "N"):?
? ? ? ? ? ? ? ? if (beg == -1):
? ? ? ? ? ? ? ? ? ? beg=pos+i+1
? ? ? ? ? ? else:?
? ? ? ? ? ? ? ? if (beg>0):
? ? ? ? ? ? ? ? ? ? print("polyN\t"+chrname+"\t"+str(beg)+"\t"+str(i+pos))
? ? ? ? ? ? ? ? ? ? beg=-1?
? ? ? ? pos = pos +len(line)
if (beg>0):
? ? print("polyN\t"+chrname+"\t"+str(beg)+"\t"+str(pos)) #如果末尾有polyN再輸出一下
測試test.fa:
>chr1
AACGNNNNACGTNNNNACGTA
ACTAGCATCGTACNNN
NNNACGTACGTACGTA
>chr2
ACTAGCATCGTACNNN
NNNACGTACGTACNNN
輸出：
#label? Chr? ? ?Beg? ? ?End
polyN? ?chr1? ? 5? ? ? ?8
polyN? ?chr1? ? 13? ? ? 16
polyN? ?chr1? ? 35? ? ? 40
polyN? ?chr2? ? 14? ? ? 19
polyN? ?chr2? ? 30? ? ? 32
這個程序?qū)懞昧耍唵问纠呀?jīng)可以獲得正確的結(jié)果，那么就可以上大基因組了。那么這時候，下一個問題就出現(xiàn)了。
4. 時間?。。?/span>
一段代碼的復(fù)雜度在各種算法課本中可以查到，不過我這里并不想講怎么計算這個復(fù)雜度。如果用上面這段代碼去跑大基因組，那么有一個問題就會出現(xiàn)，程序需要跑多久？
最簡單的方法就是，跑一遍不就知道多久了？所以一般比如要跑10G的數(shù)據(jù)，可以先拿10M或者100M的數(shù)據(jù)跑一遍，最后耗時按照比例增加就行了。
time 你的代碼.py?
耗時需要看自己的承受能力。比如10min，1h，1day，1week等等，按照自己接受能力去優(yōu)化代碼。如果已經(jīng)估計出來的時間可以接受，并且以后這段代碼不會再用，那么就沒有優(yōu)化的必要。比如耗時需要1天，但是我可能通過優(yōu)化可以變成半天，然鵝我可以去聚聚餐，打打球，明天這時候直接收結(jié)果就行了，為啥要優(yōu)化？
不過，如果估計出來的1week或者更長，自己肯定無法接受，那么就要想辦法進(jìn)行優(yōu)化。
我采用的優(yōu)化策略一般有以下2種：
用更多的資源
?優(yōu)化算法
包括原來用一臺電腦跑一周，我開7臺電腦跑一天就好了啊~或者是原來是單核的程序，可用的cpu有20核，我分上18個文件，同時跑18個不就更快；如果資源需要共享啊，交互啊什么的，可以寫多線程（不過線程鎖、線程間的共享似乎是個坑）
第一種就不說了，第二種深有體會。
4.1 算法的優(yōu)化
算法優(yōu)化很容易就可以有量級級別的提升。
故事開始：
我記得最清楚的是，我剛進(jìn)實驗室時，師兄F需要做1000w個數(shù)無重復(fù)依次抽取100w個數(shù)，他rand(1000w)后，抽一個就和抽出來的判斷一下是否抽重了，如果重了重新抽，于是這段代碼寫出來越來越慢直至卡死；之后我生成1000w的數(shù)組，之后從里面抽一個刪一個（當(dāng)時想的是從口袋中抽球），估計4-6小時可以跑完；最后，我?guī)熜諱通過數(shù)組元素交換的方法，寫的程序3秒跑完。這次對我印象太深刻了，好的算法真的可以帶來量級上的提升。
我現(xiàn)在一般用到的算法優(yōu)化策略有以下幾種：
避免循環(huán)嵌套。
避免不必要的資源浪費(fèi)。比如多余的變量，多余的數(shù)組和鏈表等。部分編程語言有較好的垃圾回收機(jī)制可以避免這種內(nèi)存占用，不過我還是習(xí)慣清空不用數(shù)組，關(guān)閉打開的文件等減少內(nèi)存占用。
數(shù)組利于按index查詢不利于遍歷搜索和增刪，鏈表利于增刪不利于遍歷，哈希鍵值查詢極快等。利用各種數(shù)據(jù)結(jié)構(gòu)的優(yōu)勢進(jìn)行代碼編寫。
引用我大學(xué)的計算機(jī)老師勇哥的話：“人想的多，計算機(jī)運(yùn)算就快”。比如1到10億的加和可以遍歷求和也可以直接用高斯的等差數(shù)列公式，很明顯后者快。
盡量把大程序拆成小塊，小的子程序，這樣也比較容易寫單元測試。
字符串查詢上，通常正則表達(dá)式要快=-=
linux比Win快=-=
linux上用shell比Perl快，文本處理Perl比Python快=-=
那么，對于上面我們寫的這個算法，我試了一下我的電腦上567M的fasta需要跑1m12s?？赡苓@個屬于我前面說的無需優(yōu)化的類型，接杯水回來就跑完了。
然鵝，如果需要跑1T的數(shù)據(jù)，就需要36.9小時，一天半。所以，按照上面的策略5，我們試試正則表達(dá)式。
4.2 上，正則表達(dá)式！
最早孟師妹問我的時候，她采用的代碼就是前面我說的（我估計是的，因為是通常的思路就是這樣的），她覺得太慢了。
于是我第一反應(yīng)就是，既然是匹配，那么上正則表達(dá)式唄。不過前提條件就是，需要內(nèi)存足夠大，序列需要放一行。如果是多行fasta正則表達(dá)式匹配，遇到3中提到的多行的問題還不夠麻煩的=-=，效率也不見得能提高多少。
所以，孟師妹修改后的代碼:
import re
print("#label\tChr\tBeg\tEnd")?
#這里直接采用孟師妹的函數(shù)，把fasta的名那行變?yōu)榕紨?shù)行，序列變?yōu)槠鏀?shù)行，首行為空
def ToSingle(fasta_file):
? ? str_my = open(fasta_file).read() #文件全部讀到內(nèi)存
? ? first_res = re.sub(r'(>.*)\n', r'@\1@', str_my) #正則表達(dá)式替換，把開頭有>也就是標(biāo)題行用兩個@包圍
? ? sec_res = first_res.replace('\n', '') #去除全部換行符（這里其實沒考慮到windows的\r）
? ? res = sec_res.replace('@', '\n') #把剛才標(biāo)記的@變成換行
? ? return res
file_arr = ToSingle(sys.argv[1]).split("\n") #用\n分割文件,生成數(shù)組
reg = re.compile(r'N+',re.I) #編譯正則表達(dá)式，多個N，忽略大小寫
chrom = ""
for index in range(1,len(file_arr)):
? ? if (file_arr[index].startswith(">")):
? ? ? ? chrom = file_arr[index].replace('>', '') #獲得染色體名
? ? else:
? ? ? ? for sub_str in reg.finditer(file_arr[index]):
? ? ? ? ? ? start, end = sub_str.span()
? ? ? ? ? ? start = start+1
? ? ? ? ? ? print(f"polyN\t{chrom}\t{start}\t{end}") #3.6新增的格式化輸出方式
這樣，我再運(yùn)行剛才567Mfasta，只用了15s，提高了4.8倍。1T數(shù)據(jù)時間可以從36.9小時減少到7.7小時。
最后思考，如果一行序列內(nèi)存放不下怎么辦？
分多行？之后再判斷切割位點(diǎn)附近序列有沒有polyN？
所以，可以看到，寫程序也有這個特點(diǎn)，有錢有資源就可以少折騰，233333。
5. 應(yīng)用性和泛性
既然我們已經(jīng)得到了一個差不多高效的程序了，還需要做什么？
當(dāng)然是整理好代碼以后還可以繼續(xù)用啊。也就是這個程序要適應(yīng)各種不同的情況，而不只是test.fa。此外，如果其他人也想用這個軟件，他們沒親自編寫程序，于是就有可能犯各種各樣的錯誤。
所以，我通常的做法就是，把上面的算法“打包”，也就是認(rèn)為“找polyA”是一個箱子，只要輸入fasta文件，他就能輸出polyN的位置。那么，我就需要把fasta文件給這個箱子而不是給個gbk的文件。
那么，就需要用到：
顯示文件用法readme，設(shè)置參數(shù)的包：argparse, sys.argv
Biopython里面完善的讀fasta文件的包: SeqIO
等。這里我就不再給出具體的代碼了。
當(dāng)然，如果想讓程序更易用，也可以加入更多的功能。比如有人需要獲得"非polyN"序列的位置，也就是把正則表達(dá)式編譯那里的
? ? reg = re.compile(r'[^N]+',re.I)
? ? #不過這里匹配非N，也許數(shù)字或者其他奇奇怪怪的也會進(jìn)來哦。
? ? #所以需要輸入數(shù)據(jù)有質(zhì)控。
等等，會有很多很多需要完善的功能，而在你完善的過程中，你的程序能力就會飛速的提升。
好啦，今天這個例子就到這里，總結(jié)一下。
這個例子主要說了我寫程序的一些思路和優(yōu)化經(jīng)驗，并展示了正則表達(dá)式。
希望大家有所收獲，歡迎與我交流。
