Linux 中 shell 脚本实现根据gff 计算每一个基因的最长转录本
001、生成基因名称的列表
awk -F "\t" '$3 == "gene" && $NF ~ /gene=/ {print $NF}' chr1.gff | sed 's/\(.*\)\(gene=[^;]\+\)\(.*\)/\2/' | sort | uniq > gene.list
002、计算最长转录本
a、脚本文件
#!/bin/bash num=0 for i in $(cat gene.list) do let num=$num+1 mRNA=$(grep "$i;" chr1.gff | awk '$3 ~ /^mRNA$/' | wc -l) exon=$(grep "$i;" chr1.gff | awk '$3 ~ /^exon$/' | wc -l) if [ $exon -eq 0 ] then echo "$i has no exon" exit fi if [ $mRNA -gt 1 ] then grep "$i;" chr1.gff | awk '$3 ~ /^mRNA$|^exon$/' | cut -f 3-5 | awk 'BEGIN{a=0}{if($1 ~ /^mRNA$/){a++}; if($1 ~ /^exon$/) {print "exon"a, $2, $3}}' | awk '{ay[$1] += ($3 - $2 + 1)} END {for(j in ay) print j, ay[j]}' | awk -v sam=$i '{if(NR == 1) {max = $2} else if($2 > max) {max = $2}} END {print sam, max}' >> result.txt else grep "$i;" chr1.gff | awk -v sam=$i '$3 ~ /^exon$/ {sum += ($5 - $4 + 1)} END {print sam, sum}' >> result.txt fi echo $num done done
b、执行程序
[root@PC1 test]# bash record.sh
c、结果文件:
[root@PC1 test]# head result.txt gene=2-CP2 1080 gene=AAP1 1684 gene=abz1 1196 gene=ACI112 1787 gene=ACI36 1445 gene=ACS2 1847 gene=AdBiL 2089 gene=ADC 2800 gene=ADH 1152 gene=Adi3 2494