搜索
查看: 4147|回复: 2

把counts矩阵转换成RPKM矩阵

[复制链接]

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2017-7-3 09:53:37 | 显示全部楼层 |阅读模式
虽然建议大家抛弃RPKM这种基因表达量的表示方式,但是因为历史原因,总是还有人需要。
overlap all exons for each gene and sum the length of each exon to obtain an estimate of the gene length
RPKMs are not the best choice to do differential analysis in gene expression.
Take a look at the following video: www.youtube.com/watch?v=TTUrtCY2k-w


所以大家用来hisat2+htseq-counts之后得到的reads的counts数量这样的矩阵,如果需要转换成RPKM矩阵

就考验编程功底啦。

我简单给一个例子哈;

[AppleScript] 纯文本查看 复制代码
wget [url=ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_mouse/CCDS.current.txt]ftp://ftp.ncbi.nlm.nih.gov/pub/C ... se/CCDS.current.txt[/url]
grep -v '^#' CCDS.current.txt | perl -alne '{/\[(.*?)\]/;$len=0;foreach(split/,/,$1){@tmp=split/-/;$len+=($tmp[1]-$tmp[0])};$h{$F[2]}=$len if $len >$h{$F[2]}} END{print "$_\t$h{$_}" foreach sort keys %h}' >mm10_ccds_length.txt



genes_len=read.table("mm10_ccds_length.txt",stringsAsFactors=F)
head(genes_len)
              V1    V2
1      -343C11.2  1139
2 00R_AC107638.2  1060
3      00R_Pgap2  1676
4  0610005C13Rik  7363
5  0610006L08Rik 34995
6  0610007P14Rik  9074

 
colnames(genes_len)<- c("GeneName","Len")

> head(exprSet)
               GSM860181_priSG-A GSM860182_SG-A GSM860183_SG-B GSM860184_lepSC
00R_AC107638.2                 0              1              0               1
0610005C13Rik                 20             22             11              27
0610006L08Rik                  0              0              0               2
0610007P14Rik               2075           1785           1269            1430
0610009B22Rik                256            376            300             386
0610009E02Rik                 22             22             16              28
 
exprSet<-exprSet[ rownames(exprSet) %in% genes_len$GeneName ,]

total_count<- colSums(exprSet)

neededGeneLength=genes_len[  match(rownames(exprSet), genes_len$GeneName) ,2] 
 
rpkm <- t(do.call( rbind,lapply(1:length(total_count),function(i){
        10^9*exprSet[,i]/neededGeneLength/total_count[i]
}) ))
head(rpkm)
rownames(rpkm)=rownames(exprSet)
colnames(rpkm)=colnames(exprSet)

write.table(rpkm,file="rpkm.txt",sep="\t",quote=F)

 
 
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb=TxDb.Mmusculus.UCSC.mm10.knownGene
gt=transcriptsBy(txdb,by="gene")
lapply(gt[1:40],function(x) max(width(x)))
library(org.Mm.eg.db)

mykeys=
columns(txdb);keytypes(txdb)
neededcols <- c("GENEID", "TXSTRAND", "TXCHROM")
select(txdb, keys=mykeys, columns=neededcols, keytype="TXNAME")





上一篇:S4
下一篇:python 如何实现自动上传文件到网页?
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

1

主题

5

帖子

117

积分

版主

Rank: 7Rank: 7Rank: 7

积分
117
发表于 2017-8-11 10:08:15 | 显示全部楼层
回复 支持 反对

使用道具 举报

0

主题

1

帖子

95

积分

注册会员

Rank: 2

积分
95
发表于 2017-9-7 23:54:27 | 显示全部楼层
Jimmy老師你好,我剛開始接觸TCGA.想請教一下如何把合並後的counts矩陣轉換成TPM矩陣?
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-21 09:36 , Processed in 0.066096 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.