RNAseq

基因长度的计算(非冗余外显子长度之和)

2021-03-12  本文已影响0人  啊吧啊吧鸭

基因的长度有以下几种定义(参考:http://www.bio-info-trainee.com/3991.html
挑选基因的最长转录本
选取多个转录本长度的平均值
非冗余外显子(EXON)长度之和
非冗余 CDS(Coding DNA Sequence) 长度之和

bash脚本(参考:https://www.jianshu.com/p/d53416962fa8

# 提取exon所在的行
awk '{if ($3=="exon") print $0}' Homo_sapiens.GRCh38.94.gtf  > exon.gtf

#去掉双引号
sed -i 's/"//g' exon.gtf

#提取需要的列到exon.txt中
awk 'BEGIN{FS="\t| |;";OFS="\t"}{print $1,$3,$4,$5,$7,$10,$25,$5-$4+1}' exon.gtf > exon.txt

#计算基因长度
awk 'BEGIN{OFS="\t"}{ s[$7] += $8 }END{for (i in s){print i,s[i]}} ' exon.txt > gene.length.txt


R脚本(参考:https://www.jianshu.com/p/4afbdc8fd8ed

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("GenomicFeatures")
library(GenomicFeatures)

# 将已有的GTF文件生成包处理的对象
txdb <- makeTxDbFromGFF("../data/F153.v20200709.gtf",format="gtf")

exons_gene <- exonsBy(txdb, by = "gene")
exons_gene_lens <- lapply(exons_gene,function(x){sum(width(reduce(x)))})
str(exons_gene_lens)
exons_gene_lens1 <- as.data.frame(exons_gene_lens)
dim(exons_gene_lens1)
exons_gene_lens1 <- t(exons_gene_lens1)

counts转TPM:https://www.jianshu.com/p/7d64b8a9fa99

上一篇 下一篇

猜你喜欢

热点阅读