1.核验病人临床药物信息是否准确（R包VS网页版Sangerbox VS XML文件一一合并）

#获取临床信息

cancer_type="TCGA-LIHC"

clinical <- GDCquery_clinic(project= cancer_type,type = "clinical")

#获取临床药物信息

query <- GDCquery(project = "TCGA-LIHC",

                  data.category = "Clinical",

                  file.type = "xml")

query_download <- GDCdownload(query)

clinical.drug <- GDCprepare_clinic(query,"drug")

#sangerbox中的Merger好的总信息TXT文件（是第一个可预览的文件）

#XML文件批量合并

XMLClinicalDataPath <- "E:/1.RStudio/tcga_array/tcga_array/Sangerbox_LIHC"

partlyXMLClinDataFN <- "partlyXMLClinData.csv"

#################

### 数据处理！！！

dir.create("tempClinXMLFiles")

filepath <- dir(path = XMLClinicalDataPath,full.names = TRUE)

### 将所有文件移动到同一个文件夹tempClinXMLFiles下

### 使用以下循环：

for(wd in filepath){

  files <- dir(path = wd,pattern = "xml$") # 查看满足条件(pattern=xml)的文件

  fromfilepath <- paste(wd,"/",files,sep = "")

  tofilepath <- paste("./tempClinXMLFiles/",files,sep = "")

  file.copy(fromfilepath,tofilepath)

}

#################

### 临床数据整合！！！

setwd("./tempClinXMLFiles")

xmlFileNames <- dir(path = "./", pattern = "xml$") # list.files()函数也可以

library(XML)

library(methods)

FromeXMLClinData <- lapply(xmlFileNames,

                          function(x){            # x = xmlFileNames[1]

                            result <- xmlParse(file = x)

                            rootnode <- xmlRoot(result)

                            xmldataframe <- xmlToDataFrame(rootnode[2])

                            return(t(xmldataframe))})

### 上面这段代码不错

FromeXMLClinData <- t(do.call(cbind, FromeXMLClinData))

### 返回上级目录，保存一下

setwd("E:/1.RStudio/tcga_array/tcga_array/Sangerbox_LIHC")

write.csv(FromeXMLClinData, file = "FromeXMLClinData.csv")

### 转换成数据框格式

FromeXMLClinData <- as.data.frame(FromeXMLClinData)

### 提取部分我们需要的临床数据：

partlyXMLClinData <- data.frame(Barcode = FromeXMLClinData$bcr_patient_barcode,

                                vital_status = FromeXMLClinData$vital_status,

                                days_to_death = FromeXMLClinData$days_to_death,

                                lastFollowupTime = FromeXMLClinData$days_to_last_followup,

                                Stage = FromeXMLClinData$stage_event,

                                Age=FromeXMLClinData$age_at_initial_pathologic_diagnosis,

                                Drug=FromeXMLClinData$drugs,

                                Grade=FromeXMLClinData$neoplasm_histologic_grade,

                                Gender=FromeXMLClinData$gender)

### 保存一下

write.csv(partlyXMLClinData,

          file = partlyXMLClinDataFN,

          row.names = F, quote = F)

2.上面的结果除非数据更新，否则结果是一致的，可以直接用sangerbox的结果，因为它非常全面

sangerbox_merger<read.table("./Merge_7156ffc3f70ec3b64dbe10e7c51a4b0a_clinical.txt",header = T,  sep ="\t",dec = ".",quote="",comment.char = "" ,na.strings = c("NA"),fill=T)

#提取感兴趣的列的临床信息

clinical_data <- data.frame(Sample_ID=sangerbox_merger$A0_Samples,

                            Age=sangerbox_merger$A17_Age,

                            Sex=sangerbox_merger$A18_Sex,

                            Height=sangerbox_merger$A19_Height,

                            OS.time=sangerbox_merger$A1_OS,

                            Weight=sangerbox_merger$A20_Weight,

                            BMI=sangerbox_merger$A21_BMI,

                            OS=sangerbox_merger$A2_Event,

                            Stage=sangerbox_merger$A6_Stage,

                            Grade=sangerbox_merger$A7_Grade,

                            Drug=sangerbox_merger$drugs.drug.drug_name)

setwd("E:/1.RStudio/tcga_array/tcga_array/LIHC_TCGA_data/cli_exp_LIHC")

save(clinical_data,file="LIHC_clinical.Rdata")

write.csv(clinical_data, file="LIHC_clinical.csv")

cli_Sor1<-clinical_data[which(clinical_data[,11]=='Sorafenib'),]

cli_Sor2<-clinical_data[which(clinical_data[,11]=='SORAFENIB'),]

cli_Sor3<-clinical_data[which(clinical_data[,11]=='Sorafenib;Sorafenib'),]

cli_Sor<-rbind(cli_Sor1,cli_Sor2)

cli_Sor<-rbind(cli_Sor,cli_Sor3)

cli_nodrug<-clinical_data[which(clinical_data[,11]==''),]

##TCGA数据处理

tcga_expr <- data.table::fread("../Data/TCGA-LIHC.htseq_fpkm.tsv.gz") %>% as.data.frame()

View(tcga_expr[1:10,])

# 更换基因名字

tcga_expr$Ensembl_ID <- sub("(ENSG\\d+)\\..*","\\1",tcga_expr$Ensembl_ID)

View(tcga_expr[1:10,])

tcga_expr$Gene <- mapIds(x = org.Hs.eg.db,

                        keys = tcga_expr$Ensembl_ID,

                        column = "SYMBOL",

                        keytype = "ENSEMBL",

                        multiVals = "first")

# 变成可以分析的表达矩阵

tcga_expr <- tcga_expr %>%

  dplyr::select(-Ensembl_ID) %>%

  filter(!is.na(Gene)) %>%

  distinct(Gene,.keep_all = T) %>%

  column_to_rownames(var = "Gene") %>%

  as.matrix()

tcga_expr <- 2^tcga_expr - 1

3.参数设置

# data.type

#下载rna-seq的counts数据

data.type = "Gene Expression Quantification"

#下载miRNA数据

data.type = "miRNA Expression Quantification"

#下载Copy Number Variation数据

data.type = "Copy Number Segment"

#workflow.type 有三种类型分别为：

#HTSeq - FPKM-UQ：FPKM上四分位数标准化值

#HTSeq - FPKM：FPKM值/表达量值

#HTSeq - Counts：原始count数

#legacy

#这个参数主要是因为TCGA数据有两个入口可以下载，

#GDC Legacy Archive 和 GDC Data Portal，

#区别主要是注释参考基因组版本不同分别是：

#GDC Legacy Archive（hg19和GDC Data Portal（hg38）。

#参数默认为FALSE，下载GDC Data Portal（hg38）。

#这里建议是，下载转录组层面的数据使用hg38，

#下载DNA层面的数据使用hg19，

#因为比如做SNP分析的时候很多数据库没有hg38版本的数据，都是hg19的。

用程序下载文件miRNA数据样本ID和患者ID的对应方法

https://cloud.tencent.com/developer/article/1802656

数据下载完以后接下来的处理

第四步：TCGAanalyze＿Preprocessing（）对数据进行预处理：使用spearman相关系数去除数据中的异常值

＃ 去除dataPrep1中的异常值，dataPrep1数据中含有肿瘤组织和正常组织的数据

＃ TCGAanalyze＿Preprocessing（object， cor．cut ＝ 0， filename ＝ NULL，width ＝ 1000， height ＝ 1000， datatype =names（assays（object））［1］）

＃ 函数功能描述：Array Array Intensity correlation （AAIC） and correlation boxplot to define outlier

dataPrep2 ＜－ TCGAanalyze＿Preprocessing（object ＝ dataPrep1，

cor．cut ＝ 0．6，datatype ＝ ＂HTSeq － Counts＂）

＃将预处理后的数据dataPrep2，写入新文件“LIHC＿dataPrep．csv”

write．csv（dataPrep2，file ＝ ＂LIHC＿dataPrep．csv＂，quote ＝ FALSE）

这里将生成一个array－array intensity correlation（AAIC）相关性热图，如下：

TCGAanalyze＿Preprocessing（）中的参数：

参数用法object来自TCGAprepare的结果cor．cut设置阈值，根据样本中各个样本之间的spearman相关系数进行过滤。默认为0filename设置生成图片文件的名称，默认为PreprocessingOutput．pngwidth生成图片的宽度 height生成图片的高度datatype描述RangedSummarizedExperiment 数据类型的字符串

第五步：TCGAtumor＿purity（）筛选肿瘤纯度大于60％的肿瘤barcodes

＃ TCGAtumor＿purity（barcodes， estimate， absolute， lump， ihc， cpe），使用来自5种方法的5个估计值作为阈值对TCGA样本进行过滤，这5个值是estimate， absolute， lump， ihc， cpe，这里设置cpe＝0．6（cpe是派生的共识度量，是将所有方法的标准含量归一化后的均值纯度水平，以使它们具有相等的均值和标准差）

＃筛选肿瘤纯度大于等于60％的样本数据

purityDATA ＜－ TCGAtumor＿purity（colnames（dataPrep1）， 0， 0， 0， 0， 0．6）

＃ filtered 为被过滤的数据， pure＿barcodes是我们要的肿瘤数据

Purity．LIHC＜－purityDATA＄pure＿barcodes

normal．LIHC＜－purityDATA＄filtered

filtered 为被过滤的数据（为正常组织的数据barcodes）， pure＿barcodes是我们要的肿瘤样本barcodes。

第六步：将肿瘤表达矩阵与正常组织表达矩阵合并，进行基因注释

＃获取肿瘤纯度大于60％的340个肿瘤组织样本＋50个正常组织样本，共计390个样本

puried＿data ＜－dataPrep2［，c（Purity．LIHC，normal．LIHC）］

第七步：进行表达矩阵基因注释

＃基因注释，需要加载“SummarizedExperiment”包，“SummarizedExperiment container”每个由数字或其他模式的类似矩阵的对象表示。行通常表示感兴趣的基因组范围和列代表样品。

＃if （！requireNamespace（＂BiocManager＂， quietly ＝ TRUE））

install．packages（＂BiocManager＂）

＃BiocManager：：install（＂SummarizedExperiment＂）    ＃没有的需要执行下载代码

library（＂SummarizedExperiment＂）

rowData（dataPrep1）   ＃传入数据dataPrep1必须为SummarizedExperiment对象

＃ DataFrame with 56512 rows and 3 columns

＃                 ensembl＿gene＿id external＿gene＿name original＿ensembl＿gene＿id

＃                     ＜character＞        ＜character＞              ＜character＞

＃ ENSG00000000003 ENSG00000000003             TSPAN6       ENSG00000000003．13

＃ ENSG00000000005 ENSG00000000005               TNMD        ENSG00000000005．5

＃ ENSG00000000419 ENSG00000000419               DPM1       ENSG00000000419．11

＃ ENSG00000000457 ENSG00000000457              SCYL3       ENSG00000000457．12

＃将结果写入文件“puried．LIHC．cancer．csv”

rownames（puried＿data）＜－rowData（dataPrep1）＄external＿gene＿name

write．csv（puried＿data，file ＝ ＂puried．LIHC．csv＂，quote ＝ FALSE）

第八步：进行表达矩阵标准化和过滤，得到用于差异分析的表达矩阵

｀TCGAanalyze＿Normalization（）｀使用EDASeq软件包标准化mRNA转录本和miRNA。

＃TCGAanalyze＿Normalization（）执行EDASeq包中的如下功能：

1． EDASeq：：newSeqExpressionSet

2． EDASeq：：withinLaneNormalization

3． EDASeq：：betweenLaneNormalization

4． EDASeq：：counts

dataNorm ＜－ TCGAanalyze＿Normalization（tabDF ＝ puried＿data，

geneInfo ＝ geneInfo，

method ＝ ＂gcContent＂）

TCGAanalyze＿Normalization中的参数：

参数用法tabDFRNAseq表达矩阵，行代表基因，列代表样本geneInfo关于geneLength和gcContent的20531个基因的矩阵，“geneInfoHT”和“geneInfo”可选。method选择标准化的方法，基于’gcContent’ 或 ’geneLength’的标准化方法可选

缺失值补全写作：

首先我们获取了表达矩阵，去除了NA值比例大于X%的行（基因）和NA值比例大于Y%的列（样本），更进一步的利用R软件包impute的impute.knn函数进行缺失值补全，设置Number of neighbors 为 K以补全数据缺失。

数据标准化部分写作：

更进一步的我们对数据进行了【分位数标准化、中位数标准化、log2转换，行Z-score转换，列Z-score转换，行和列Z-score转换

＃将标准化后的数据再过滤，去除掉表达量较低（count较低）的基因，得到最终的数据

dataFilt ＜－ TCGAanalyze＿Filtering（tabDF ＝ dataNorm，

method ＝ ＂quantile＂，

qnt．cut ＝  0．25）

str（dataFilt）

＃num ［1：13083， 1：340］ 274 2432 60347 1012 1947 ．．．

＃－ attr（＊， ＂dimnames＂）＝List of 2

＃ ．．＄ ： chr ［1：13083］ ＂A1BG＂ ＂A1CF＂ ＂A2M＂ ＂A4GALT＂ ．．．

＃ ．．＄ ： chr ［1：390］ ＂TCGA－DD－AAD5－01A－11R－A41C－07＂ ＂TCGA－DD－A4NO－01A－11R－A28V－07＂ ＂TCGA－EP－A2KA－01A－11R－A180－07＂ ＂TCGA－DD－AACP－01A－11R－A41C－07＂ ．．．

TCGAanalyze＿Filtering（）中的参数：

参数用法tabDF数据框或者矩阵，行代表基因，列代表来自TCGA的样本method用于过滤较低count数的基因的方法，有’quantile’， ’varFilter’， ’filter1’， ’filter2’qnt．cut选择均值作为过滤的阈值

最后将过滤后的数据写入文件“TCGA＿LIHC＿final．csv”，就得到我们用于后续差异分析的表达文件：

write．csv（dataFilt，file ＝ ＂TCGA＿LIHC＿final．csv＂，quote ＝ FALSE）

＃保留的是390个样本（前340肿瘤，后50正常组织）