分享一個(gè)基因組數(shù)據(jù)篩選過(guò)程中遇到的簡(jiǎn)單案例

我有一個(gè)gtf文件，格式如下圖1所示；我想將每一行中g(shù)ene_id部分都篩選出來(lái)（如圖1中紅圈）。
但是首要的問(wèn)題是下載的gtf文件，有的行中可能沒(méi)有g(shù)ene_id，因此我想做的首先就是判斷每一行中是否都有g(shù)ene_id，如果有，則判斷為True，如果沒(méi)有，則判斷為False，這一步可寫一個(gè)python腳本實(shí)現(xiàn)；

image1.png

判斷每一行中是否“gene_id” 這一字符串

#!/usr/bin/env python
inputfile = "GCF_000001735.4_TAIR10.1_genomic.gtf"
outputfile = "yangjincheng.txt"
output = open(outputfile,"w")
line_count = 1
for line in open(inputfile):
        output.write(str(line_count) + "\t" + str("gene_id" in line) + "\n")      ###最重要的是 A in line 這個(gè)運(yùn)算符，具有判斷功能
        line_count += 1
output.close(）

得到的結(jié)果如圖2所示：第一列是行名，第二列是判斷結(jié)果。因此我們就需要把第二列中顯示是True的行篩選出來(lái)；

image2.png

篩選判斷結(jié)果是True的行

sed "/True/p" yangjincheng.txt -n | awk '{print $1}' > yangjincheng1.txt         ### 用三劍客中的sed和awk就能篩掉判斷結(jié)果為False的行名

得到結(jié)果如圖3所示；得到了含有判斷結(jié)果為True的行名，接下來(lái)就需要將對(duì)應(yīng)行名的行從gtf文件中提取出來(lái)，可以用sed命令實(shí)現(xiàn)，但為了批量處理，應(yīng)首先寫一個(gè)python腳本；

image3.png

利用python腳本寫出sed命令的標(biāo)準(zhǔn)輸出集合

#!/usr/bin/env python
inputfile = "yangjincheng1.txt"
outputfile = "jincheng.sh"
output = open(outputfile,"w")
for line in open(inputfile):
        output.write('sed ' + '"' + line.rstrip() + 'p' + '" ' + 'GCF_000001735.4_TAIR10.1_genomic.gtf ' + '-n' + '\n')
output.close()

利用該腳本生成了一個(gè)sed命令的集合，如圖4所示；其實(shí)是生成了一個(gè)shell腳本，能實(shí)現(xiàn)數(shù)據(jù)的批量化處理

image4.png

運(yùn)行新生成的shell腳本

sh jincheng.sh > GCF.gtf &         ###運(yùn)行速度比較慢，可能有更好的解決辦法

這樣就把最初的gtf文件中，行中含有“gene_id”的數(shù)據(jù)都篩選了出來(lái)，下一步將利用一個(gè)python腳本將GCF.gtf中的gene_id部分篩選出來(lái)。

篩選出所有的gene_id

#!/usr/bin/env python
inputfile = "GCF.gtf"
outputfile = "chenyu.txt"
output = open(outputfile,"w")
output.write("gene_id" + "\n")
for line in open(inputfile):
        line1 = line[(line.index("gene_id") + 9):-1]
        gene_id = line1[:line1.index('"')]
        output.write(gene_id + "\n")
output.close()

運(yùn)行最后這個(gè)腳本就能得到想要的結(jié)果