搜索
查看: 6024|回复: 3

把counts矩阵转换成RPKM矩阵

[复制链接]

633

主题

1189

帖子

4054

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4054
发表于 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矩阵

就考验编程功底啦。

我简单给一个例子哈;

[mw_shl_code=applescript,true]
wget ftp://ftp.ncbi.nlm.nih.gov/pub/C ... se/CCDS.current.txt
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
}) ))
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")
[/mw_shl_code]





上一篇: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矩陣?
回复 支持 反对

使用道具 举报

0

主题

1

帖子

118

积分

注册会员

Rank: 2

积分
118
发表于 2020-1-5 15:57:50 | 显示全部楼层
计算RPKM那一步是不是有个错误,就是在那一行total_count应该写成total_count[i]
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2020-4-5 03:43 , Processed in 0.024532 second(s), 24 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.