bash脚本计算转录本长度

转录本长度，在计算基因表达量TPM等需要，本脚本从gtf文件提取生成转录本长度文件。

1. 找出所有的转录本

# 找出所有的转录本
cat hg38.refGene.gtf|awk -F '\t' '{if($3=="transcript")  print $0}'|head

2. 找出某一转录本的所有外显子

# 某一转录本的所有外显子
cat hg38.refGene.gtf|grep "\"NM_001476\""|awk '{if($3=="exon")  print $0}' # 匹配"NM_001476"

3. 计算某一转录本所有外显子之和

# 计算某一转录本的长度（所有exon之和）
cat hg38.refGene.gtf|grep "\"NM_001476\""|awk '{if($3=="exon")  print $0}'|awk '{print $5,$4,$5-$4+1}'|awk '{sum+=$3};END {print "NM_001476\t" sum}'

4. 找出某一基因的所有转录本

# 找出某一基因的所有转录本，如BRAF
cat hg38.refGene.gtf|awk -F '\t' '{if($3=="transcript")  print $0}'|grep BRAF

5. 列出所有转录本的编码区长度（外显子之和）

## 列出转录本所有exon之和
##分两步
cat hg38.refGene.gtf|awk -F '\t' '{if($3=="transcript")  print $9}'|awk '{print $4}' > transtript_id.txt

for id in `cat transtript_id.txt|head -2`; do cat hg38.refGene.gtf|grep $id|awk '{if($3=="exon") print $0}'|awk '{print $5,$4,$5-$4+1}'|awk -v name=${id:1:-2} '{sum+=$3};END {print name, "\t", sum}';done

# 注意-v参数 ${id:1:-2} 去掉前面一个和后面两个字符

## 一行代码搞定
for id in `cat hg38.refGene.gtf|awk -F '\t' '{if($3=="transcript")  print $9}'|awk '{print $4}'`; do cat hg38.refGene.gtf|grep $id|awk '{if($3=="exon") print $0}'|awk '{print $5,$4,$5-$4+1}'|awk -v name=${id:1:-2} '{sum+=$3};END {print name, "\t", sum}';done  >cds_len.txt

6. 某一转录本mRNA长度

cat hg38.refGene.gtf|grep "\"NM_001476\""|awk '{if($3=="exon"||$3=="5UTR"||$3=="3UTR") print $0}'|awk '{print $5,$4,$5-$4+1}'|awk '{sum+=$3};END {print "NM_001476\t" sum}'

7. 列出所有转录本mRNA长度

for id in `cat hg38.refGene.gtf|awk -F '\t' '{if($3=="transcript")  print $9}'|awk '{print $4}'`; do cat hg38.refGene.gtf|grep $id|awk '{if($3=="exon" || $3=="5UTR" || $3=="3UTR") print $0}'|awk '{print $5,$4,$5-$4+1}'|awk -v name=${id:1:-2} '{sum+=$3};END {print name, "\t", sum}';done

注：某一转录本加工成成熟的mRNA,长度应该为 5UTR+CDS+3UTR