分享一个基因组数据筛选过程中遇到的简单案例
我有一个gtf文件，格式如下图1所示；我想将每一行中gene_id部分都筛选出来（如图1中红圈）。
但是首要的问题是下载的gtf文件，有的行中可能没有gene_id，因此我想做的首先就是判断每一行中是否都有gene_id，如果有，则判断为True，如果没有，则判断为False，这一步可写一个python脚本实现；
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 这个运算符，具有判断功能
line_count += 1
output.close(）
得到的结果如图2所示：第一列是行名，第二列是判断结果。因此我们就需要把第二列中显示是True的行筛选出来；
image2.png
筛选判断结果是True的行
sed "/True/p" yangjincheng.txt -n | awk '{print $1}' > yangjincheng1.txt ### 用三剑客中的sed和awk就能筛掉判断结果为False的行名
得到结果如图3所示；得到了含有判断结果为True的行名，接下来就需要将对应行名的行从gtf文件中提取出来，可以用sed命令实现，但为了批量处理，应首先写一个python脚本；
image3.png
利用python脚本写出sed命令的标准输出集合
#!/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()
利用该脚本生成了一个sed命令的集合，如图4所示；其实是生成了一个shell脚本，能实现数据的批量化处理
image4.png
运行新生成的shell脚本
sh jincheng.sh > GCF.gtf & ###运行速度比较慢，可能有更好的解决办法
这样就把最初的gtf文件中，行中含有“gene_id”的数据都筛选了出来，下一步将利用一个python脚本将GCF.gtf中的gene_id部分筛选出来。
筛选出所有的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()
运行最后这个脚本就能得到想要的结果
