TCGA数据挖掘（四）：表达差异分析（4）

setwd("E:\\BioInfoCloud\\TCGABiolinks\\case_study1")
# 加载相应的包，可能会需要其他包，提示错误就安装缺少的包。
# 因为每个人已经安装的包不一样。
library(TCGAbiolinks)
library(plyr)
library(limma)
library(biomaRt)
library(SummarizedExperiment)
# 请求数据。
query <- GDCquery(project = "TCGA-BRCA",
                  data.category = "Transcriptome Profiling",
                  data.type = "Gene Expression Quantification", 
                  workflow.type = "HTSeq - Counts")
# 从query中获取结果表，它可以选择带有cols参数的列，并使用rows参数返回若干行。
# 1222个barcode                 
samplesDown <- getResults(query,cols=c("cases"))
# samplesDown中筛选出TP(主要实体瘤)样本的barcode
# 1102个barcode
dataSmTP <- TCGAquery_SampleTypes(barcode = samplesDown,
                                  typesample = "TP")
# 113个barcode
dataSmNT <- TCGAquery_SampleTypes(barcode = samplesDown,
                                  typesample = "NT")

# 选择100个正常组织和100个肿瘤组织样本作为研究对象
dataSmTP_short <- dataSmTP[1:100]
dataSmNT_short <- dataSmNT[1:100]
# 根据前面的筛选，再次请求数据
queryDown <- GDCquery(project = "TCGA-BRCA", 
                      data.category = "Transcriptome Profiling",
                      data.type = "Gene Expression Quantification", 
                      workflow.type = "HTSeq - Counts", 
                      barcode = c(dataSmTP_short, dataSmNT_short))
#读取下载的数据并将其准备到R对象中，在工作目录生成brca_case1.rda文件，
#同时还产生Human_genes__GRCh38_p12_.rda文件
GDCdownload(query = queryDown)
dataPrep1 <- GDCprepare(query = queryDown, save = TRUE, save.filename = "brca_case1.rda")

# 去除dataPrep1中的异常值，dataPrep数据中含有肿瘤组织和正常组织的数据
dataPrep <- TCGAanalyze_Preprocessing(object = dataPrep1, 
                                      cor.cut = 0.6,
                                      datatype = "HTSeq - Counts")
#使用EDASeq进行标准化
dataNorm <- TCGAanalyze_Normalization(tabDF = dataPrep,
                                      geneInfo = geneInfoHT,
                                      method = "gcContent")
#根据分位数一处count较低的基因。
dataFilt <- TCGAanalyze_Filtering(tabDF = dataNorm,
                                  method = "quantile", 
                                    qnt.cut =  0.25)
# 将计数数据转换为log2计数/百万(LogCPM)，估计均值-方差关系并使用此关系计算适当的观测级别权重。
# 然后，数据即可用于线性建模。
v.dataFilt<-voom(dataFilt)
#in order to have 2 batches with multiple samples and avoid batches with one sample
#the user need to call first the get_IDs function on the top of this script if this has not been done already
c1<-which(get_IDs(dataPrep)$tss=="E9")
c2<-which(get_IDs(dataPrep)$tss=="E2")
#taking log transformed data for exploration of batch effects
TCGAbatch_Correction(tabDF = v.dataFilt$E[,c(c1,c2)], batch.factor="TSS", adjustment=NULL)
#if you get error in plot.new():figure margins too large - please use:
#par(mar=c(1,1,1,1))
#and then rerun the command above

#### After exploration of the batches we can continue to work on the original gene matrix for DEA ###
###############Tumor purity filtering###########
###vector containing all TCGA barcodes that hhave 60% tumor purity or more
Purity.BRCA<-TCGAtumor_purity(colnames(dataPrep1), 0, 0, 0, 0, 0.6)$pure_barcodes
################DEA with Molecular subtypes############
####Molecular subtypes:
####Filtering data so all samples have a pam50 subtype for BRCA:

#diff contains TCGA samples that have an available molecular subtype
###Also Applying Tumor purity filtering by making a set difference
#documentation for TCGA_MolecularSubtype is available
##########Important########
####in this example no change happens because all the TCGA barcodes
####have available subtype info
###setdiff() is here used to remove samples with no subtype info from the TCGA BRCA samples that satisfy
###tumor purity criteria
diff<-setdiff(Purity.BRCA, TCGA_MolecularSubtype(colnames(dataPrep[,dataSmTP_short]))$filtered)
dataFilt.brca.cancer<-dataPrep[,diff]
dataFilt.brca.normal<-dataPrep[,dataSmNT_short]
dataFilt.brca<-cbind(dataFilt.brca.normal, dataFilt.brca.cancer)
mol_subtypes<-c(rep("Normal", 100), TCGA_MolecularSubtype(colnames(dataFilt.brca.cancer))$subtypes$subtype)
# All barcodes have available molecular subtype info
mol_subtypes<-make.names(mol_subtypes)
# to convert ENSEMBL id to gene name
rownames(dataFilt.brca)<-rowData(dataPrep1)$external_gene_name
# we redo here the normalization and filtering steps 
dataNorm.brca <- TCGAanalyze_Normalization(tabDF = dataFilt.brca,
                                           geneInfo = geneInfo,
                                           method = "gcContent")
# 将标准化后的数据再过滤，得到最终的数据
dataFilt.brca.final <- TCGAanalyze_Filtering(tabDF = dataNorm.brca,
                                       method = "quantile", 
                                       qnt.cut =  0.25)

DEG.brca.TSS <- TCGAanalyze_DEA(MAT=dataFilt.brca.final,
                            pipeline="limma",
                            batch.factors = c("TSS"),
                            Cond1type = "Normal",
                            Cond2type = "Tumor",
                            fdr.cut = 0.01 ,
                            logFC.cut = 1,
                            voom = TRUE,
                            Condtypes = mol_subtypes,
                            contrast.formula = "Mycontrast = (BRCA.LumA+BRCA.LumB)/2 -Normal")
#in this DEA we use Normal as a reference, thus genes with LogFC > 1 are up regulated in the subtypes
#with respect to the normal samples and viceversa for the LogFC < -1.

write.csv(DEG.brca.TSS, "DEGsBRCA_limma_091018.csv", quote = FALSE)
#3241 genes identified

#to check how many with logFC > 5 (for plotting purposes)
DEG.brca.filt.TSS<-DEG.brca.TSS$Mycontrast[which(abs(DEG.brca.TSS$Mycontrast$logFC) >= 6), ]
#47 genes identified