2022-09-18

手机游戏开发者 · 2024-9-5 05:53:45

分享一个基因组数据筛选过程中遇到的简朴案例

我有一个gtf文件，格式如下图1所示；我想将每一行中gene_id部分都筛选出来（如图1中红圈）。
但是主要的问题是下载的gtf文件，有的行中大概没有gene_id，因此我想做的起首就是判定每一行中是否都有gene_id，如果有，则判定为True，如果没有，则判定为False，这一步可写一个python脚本实现；

判定每一行中是否“gene_id” 这一字符串

#!/usr/bin/env pythoninputfile = "GCF_000001735.4_TAIR10.1_genomic.gtf"outputfile = "yangjincheng.txt"output = open(outputfile,"w")line_count = 1for line in open(inputfile):       output.write(str(line_count) + "\t" + str("gene_id" in line) + "\n")    ###最紧张的是 A in line 这个运算符，具有判定功能       line_count += 1output.close(）得到的结果如图2所示：第一列是行名，第二列是判定结果。因此我们就必要把第二列中表现是True的行筛选出来；

筛选判定结果是True的行

sed "/True/p" yangjincheng.txt -n | awk '{print $1}' > yangjincheng1.txt       ### 用三剑客中的sed和awk就能筛掉判定结果为False的行名得到结果如图3所示；得到了含有判定结果为True的行名，接下来就必要将对应行名的行从gtf文件中提取出来，可以用sed下令实现，但为了批量处理处罚，应起首写一个python脚本；

使用python脚本写出sed下令的尺度输出集合

#!/usr/bin/env pythoninputfile = "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脚本，能实现数据的批量化处理处罚

运行新天生的shell脚本

sh jincheng.sh > GCF.gtf &       ###运行速率比较慢，大概有更好的解决办法如许就把最初的gtf文件中，行中含有“gene_id”的数据都筛选了出来，下一步将使用一个python脚本将GCF.gtf中的gene_id部分筛选出来。
筛选出全部的gene_id

#!/usr/bin/env pythoninputfile = "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()运行末了这个脚本就能得到想要的结果

2022-09-18

所属分类: 问答交流

新帖推荐: 30日

推荐作品