2021-07-04史上最暴力提取轉(zhuǎn)錄本和基因長(zhǎng)度(基于shell命令行)

在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

最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成,瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明:文章內(nèi)容(如有圖片或視頻亦包括在內(nèi))由作者上傳并發(fā)布,文章內(nèi)容僅代表作者本人觀點(diǎn),簡(jiǎn)書系信息發(fā)布平臺(tái),僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容