生存分析代码的原始版本见：两种方法批量生存分析，这个是R包GDCRNATools给出的简化版。大段代码变成几个函数搞定。
0.R包和数据准备
if(!require(GDCRNATools))BiocManager::install("GDCRNATools")
library(GDCRNATools)
# mRNA 和miRNA的表达矩阵
data(rnaCounts);dim(rnaCounts)
## [1] 1000 45
rnaCounts[1:3,1:3]
## TCGA-3X-AAV9-01 TCGA-3X-AAVA-01
## ENSG00000003989 1520 960
## ENSG00000004799 7659 957
## ENSG00000005812 2246 1698
## TCGA-3X-AAVB-01
## ENSG00000003989 2177
## ENSG00000004799 2295
## ENSG00000005812 2454
data(mirCounts);dim(mirCounts)
## [1] 2588 45
mirCounts[1:3,1:3]
## TCGA-3X-AAV9-01 TCGA-3X-AAVA-01
## hsa-let-7a-5p 165141 132094
## hsa-let-7a-3p 204 169
## hsa-let-7a-2-3p 30 26
## TCGA-3X-AAVB-01
## hsa-let-7a-5p 210259
## hsa-let-7a-3p 298
## hsa-let-7a-2-3p 50
#临床信息
metaMatrix.RNA <- gdcParseMetadata(project.id = 'TCGA-CHOL',
data.type = 'RNAseq',
write.meta = FALSE)
metaMatrix.RNA <- gdcFilterDuplicate(metaMatrix.RNA)
metaMatrix.RNA <- gdcFilterSampleType(metaMatrix.RNA)
metaMatrix.RNA[1:4,1:4]
## file_name
## TCGA-3X-AAV9-01A 725eaa94-5221-4c22-bced-0c36c10c2c3b.htseq.counts.gz
## TCGA-3X-AAVA-01A b6a2c03a-c8ad-41e9-8a19-8f5ac53cae9f.htseq.counts.gz
## TCGA-3X-AAVB-01A c2765336-c804-4fd2-b45a-e75af2a91954.htseq.counts.gz
## TCGA-3X-AAVC-01A 8b20cba8-9fd5-4d56-bd02-c6f4a62767e8.htseq.counts.gz
## file_id
## TCGA-3X-AAV9-01A 85bc7f81-51fb-4446-b12d-8741eef4acee
## TCGA-3X-AAVA-01A 42b8d463-6209-4ea0-bb01-8023a1302fa0
## TCGA-3X-AAVB-01A 6e2031e9-df75-48df-b094-8dc6be89bf8b
## TCGA-3X-AAVC-01A 19e8fd21-f6c8-49b0-aa76-109eef46c2e9
## patient sample
## TCGA-3X-AAV9-01A TCGA-3X-AAV9 TCGA-3X-AAV9-01
## TCGA-3X-AAVA-01A TCGA-3X-AAVA TCGA-3X-AAVA-01
## TCGA-3X-AAVB-01A TCGA-3X-AAVB TCGA-3X-AAVB-01
## TCGA-3X-AAVC-01A TCGA-3X-AAVC TCGA-3X-AAVC-01
rnaExpr <- gdcVoomNormalization(counts = rnaCounts, filter = FALSE)
mirExpr <- gdcVoomNormalization(counts = mirCounts, filter = FALSE)
1.差异分析
table(metaMatrix.RNA$sample_type)
##
## PrimaryTumor SolidTissueNormal
## 36 9
DEGAll <- gdcDEAnalysis(counts = rnaCounts,
group = metaMatrix.RNA$sample_type,
comparison = 'PrimaryTumor-SolidTissueNormal',
method = 'limma');dim(DEGAll)
## [1] 565 8
head(DEGAll)
## symbol group logFC AveExpr
## ENSG00000143257 NR1I3 protein_coding -6.916825 7.023129
## ENSG00000205707 ETFRF1 protein_coding -2.492182 9.515997
## ENSG00000134532 SOX5 protein_coding -4.871118 6.228227
## ENSG00000141338 ABCA8 protein_coding -5.653794 7.520581
## ENSG00000066583 ISOC1 protein_coding -2.370131 10.466194
## ENSG00000164188 RANBP3L protein_coding -5.624376 4.356284
## t PValue FDR B
## ENSG00000143257 -17.29086 4.244355e-22 2.419282e-19 40.04288
## ENSG00000205707 -16.06753 8.353256e-21 2.380678e-18 37.19751
## ENSG00000134532 -15.03589 1.168746e-19 2.220617e-17 34.49828
## ENSG00000141338 -14.86069 1.851519e-19 2.638414e-17 34.11581
## ENSG00000066583 -14.56532 4.053959e-19 4.621513e-17 33.35640
## ENSG00000164188 -14.22477 1.013592e-18 9.629125e-17 32.25659
可以获取全部差异基因，也可以单独获取mRNA和lncRNA的差异分析结果
deALL <- gdcDEReport(deg = DEGAll, gene.type = 'all');dim(deALL)
## [1] 283 8
deLNC <- gdcDEReport(deg = DEGAll, gene.type = 'long_non_coding');dim(deLNC)
## [1] 47 8
dePC <- gdcDEReport(deg = DEGAll, gene.type = 'protein_coding');dim(dePC)
## [1] 222 8
2.任意两个基因的相关性图
gdcCorPlot(gene1 = 'ENSG00000003989',
gene2 = 'ENSG00000004799',
rna.expr = rnaExpr,
metadata = metaMatrix.RNA)
3.生存分析
支持两种方法：CoxPH和KM，基于survival包，函数是gdcSurvivalAnalysis()。
CoxPH analysis
####### CoxPH analysis #######
survOutput <- gdcSurvivalAnalysis(gene = rownames(deALL),
method = 'coxph',
rna.expr = rnaExpr,
metadata = metaMatrix.RNA)
table(survOutput$pValue<0.05)
KM analysis
####### KM analysis #######
survOutput <- gdcSurvivalAnalysis(gene = rownames(deALL),
method = 'KM',
rna.expr = rnaExpr,
metadata = metaMatrix.RNA,
sep = 'median')
table(as.numeric(as.character(survOutput$pValue))<0.05)
KM plot
gdcKMPlot(gene = 'ENSG00000003989',
rna.expr = rnaExpr,
metadata = metaMatrix.RNA,
sep = 'median')
