在count矩陣轉(zhuǎn)化成TPM或者fpkm矩陣時(shí)會(huì)用到基因或者轉(zhuǎn)錄本長(zhǎng)度進(jìn)行均一化處理,那基因和轉(zhuǎn)錄本長(zhǎng)度具體怎么計(jì)算那,寫腳本有困難?龍哥,教你史上最簡(jiǎn)單,傻瓜操作!
1.基于gtf或者gff提取每個(gè)轉(zhuǎn)錄本所有exons
awk '{if ($3=="exon") print $0}' Oryza_sativa_Japonica.MSU_v7.0.gff3 > exon1.gtf

2.提取exon長(zhǎng)度及對(duì)應(yīng)轉(zhuǎn)錄本;獲取每個(gè)轉(zhuǎn)錄本所包含每個(gè)exon的長(zhǎng)度。
cut -f 4,5,9 exon1.gtf|awk '{print $3"\t"$2-$1}' >?exon2.gtf

3.計(jì)算每個(gè)轉(zhuǎn)錄本總長(zhǎng)度
awk 'BEGIN{OFS="\t"}{ s[$1] += $2 }END{for (i in s){print i,s[i]}} ' exon2.gtf >?Transcrip.length

4.生成gene - 轉(zhuǎn)錄本 -轉(zhuǎn)錄本長(zhǎng)度表格
awk -F '.' '{print $1"\t"$0}' Transcrip.length >?Transcrip.length1

5.提取gene長(zhǎng)度,gene長(zhǎng)度等于gene的最長(zhǎng)轉(zhuǎn)錄本長(zhǎng)度!
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