2021-07-04史上最暴力提取转录本和基因长度(基于shell命令行)

在count矩阵转化成TPM或者fpkm矩阵时会用到基因或者转录本长度进行均一化处理,那基因和转录本长度具体怎么计算那,写脚本有困难?龙哥,教你史上最简单,傻瓜操作!

1.基于gtf或者gff提取每个转录本所有exons

awk '{if ($3=="exon") print $0}' Oryza_sativa_Japonica.MSU_v7.0.gff3 > exon1.gtf


2.提取exon长度及对应转录本;获取每个转录本所包含每个exon的长度。

cut -f 4,5,9 exon1.gtf|awk '{print $3"\t"$2-$1}' > exon2.gtf


3.计算每个转录本总长度

awk 'BEGIN{OFS="\t"}{ s[$1] += $2 }END{for (i in s){print i,s[i]}} ' exon2.gtf > Transcrip.length


4.生成gene - 转录本 -转录本长度表格

awk -F '.' '{print $1"\t"$0}' Transcrip.length > Transcrip.length1


5.提取gene长度,gene长度等于gene的最长转录本长度!

perl -e '{while(<>){chomp;@A=split/\t/;$h{$A[0]} =$A[2] if $A[2] >$h{$A[0]}}} END{print "$_\t$h{$_}\n" foreach sort keys %h}' Transcrip.length1 > gene.length

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容