探索并获取clinical表格和miRNA表达矩阵。
参考学习资料：简书作者小洁忘了怎么分身TCGA2.GDC数据整理
1.xml文件探索
使用R包XML：
#install.packages("XML")
library(XML)
result <- xmlParse("~/Downloads/clinical/014afcf0-da24-4c9b-9017-c745f04059ae/nationwidechildrens.org_clinical.TCGA-ZC-AAAF.xml")
rootnode <- xmlRoot(result)
rootsize <- xmlSize(rootnode)
print(rootnode[1])
print(rootnode[2])
初步探索发现，xml共有两个节点，第二个节点中储存着病人的信息。
> print(rootnode[1])
$admin
...
> print(rootnode[2])
$patient
...
进一步探索使用函数xmlToDataFrame获取临床信息相关内容
xmldataframe <- xmlToDataFrame(rootnode[2])
t(xmlToDataFrame(rootnode[2]))
> t(xmlToDataFrame(rootnode[2]))
[,1]
additional_studies ""
tissue_source_site "ZC"
patient_id "AAAF"
bcr_patient_barcode "TCGA-ZC-AAAF"
bcr_patient_uuid "CD458290-F27B-479C-9B78-0AF23E7F63D9"
informed_consent_verified "YES"
icd_o_3_site "C38.1"
icd_o_3_histology "8582/3"
icd_10 "C38.1"
tissue_prospective_collection_indicator "YES"
tissue_retrospective_collection_indicator "NO"
days_to_birth "-27099"
gender "FEMALE"
height "165"
weight "68"
race_list "WHITE"
ethnicity "NOT HISPANIC OR LATINO"
other_dx "No"
history_of_neoadjuvant_treatment "No"
person_neoplasm_cancer_status "TUMOR FREE"
vital_status "Alive"
days_to_last_followup "874"
days_to_death ""
radiation_therapy "NO"
postoperative_rx_tx "NO"
post_op_ablation_embolization_tx "NO"
stage_event "IIb"
primary_pathology "Anterior MediastinumThymoma; Type ABMedian Sternectomy0742011YES"
new_tumor_events "NO"
day_of_form_completion "21"
month_of_form_completion "5"
year_of_form_completion "2014"
follow_ups "TCGA-ZC-AAAF-F685785E4486F6-267F-4171-861A-1F7184161FD3NONONOTUMOR FREEAlive1165NO3122014TCGA-ZC-AAAF-F70469563EAF41-4576-4D4C-95FB-4722FFF6B0C1NONONOTUMOR FREEAlive1165NO1822015"
drugs ""
radiations ""
病例相关信息记载非常详细，这个病人是个白人女性，无肿瘤，疾病编码是C38.1，已经随访了874天，具体信息见上述信息。
其中一个病人的信息是这样获取的，总共下载了116个病例信息，需要循环代码来读取。
学习小洁老师的循环代码
xmls = dir("~/Downloads/clinical/",pattern = "*.xml$",recursive = T)
td = function(x){
result <- xmlParse(file.path("~/Downloads/clinical/",x))
rootnode <- xmlRoot(result)
xmldataframe <- xmlToDataFrame(rootnode[2])
return(t(xmldataframe))
}
cl = lapply(xmls,td)
cl_df <- t(do.call(cbind,cl))
cl_df[1:3,1:3]
length(cl)
dim(cl[[1]])
dim(cl[[2]])
在这我到了cl_df <- t(do.call(cbind,cl))这一步的时候出错了，卡了很久，我尝试了很多种办法一直搞不明白哪里出了问题。
后来想到一种可能就是在前面我在第一个环节选择的时候在project这个选项选了2个，估计是这个地方出错，导致我得到的临床信息的记录是不一致的，所以出现了合并失败的结果。
出错原因查找
这里分析一下为什么会出现这个错误。又该怎么解决？
首先在TCGA数据库里面关于心脏相关的样本组织很少，一个是因为心脏的样本比较少，我当时下意识的选了2个项目，估计是想样本量多一点，但是不应该在方法还不熟练的情况下这样选择
解决方案有2种，一个是从新来一遍直选一个项目重头开始下载3个数据文件在一个一个合并。另一个是写一个算法来把所有四个项目下的样本都提取出来，比如只提取样本信息的某几项共有的指标，比如年龄性别，疾病，随访情况等，这样的话合并起来应该会容易一点，但是由于我的水平还达不到，暂时不能完成。
另一个就是既然我遇到了这种情况，那么其他人应该也遇到过，估计别人已经写好了，我只要找到这个方法即可。
2.探索miRNA信息
> x = read.table(file = "~/Downloads/mirna/032fd9b2-2400-4bce-b105-1f56def3257d/527618eb-b661-4f85-b586-4f2042f4f509.mirbase21.mirnas.quantification.txt",header = T,sep = "\t")
> head(x)
miRNA_ID read_count reads_per_million_miRNA_mapped cross.mapped
1 hsa-let-7a-1 108784 12676.0222 N
2 hsa-let-7a-2 109625 12774.0195 N
3 hsa-let-7a-3 109658 12777.8648 N
4 hsa-let-7b 148125 17260.2201 N
5 hsa-let-7c 20287 2363.9364 N
6 hsa-let-7d 3243 377.8896 N
> tail(x)
miRNA_ID read_count reads_per_million_miRNA_mapped cross.mapped
1876 hsa-mir-95 29 3.379216 N
1877 hsa-mir-9500 0 0.000000 N
1878 hsa-mir-96 160 18.643951 N
1879 hsa-mir-98 428 49.872569 N
1880 hsa-mir-99a 5991 698.099436 Y
1881 hsa-mir-99b 167587 19528.023723 N
可以看到，数据信息记载的第一列是miRNA_ID,第二列是raw_counts,第三列是RPM，至于第四列是什么没搞明白，但是目前重要的是我们需要前面两列的内容。
同样还是套用小洁老师总结的循环代码来汇总miRNA信息
mis = dir("~/Downloads/mirna/",pattern = "*mirnas.quantification.txt$",recursive = T)
ex = function(x){
result <- read.table(file.path("~/Downloads/mirna/",x),sep = "\t",header = T)[,1:2]
return(result)
}
mi = lapply(mis,ex)
mi_df <- t(do.call(cbind,mi))
mi_df[1:4,1:4]
得到表达矩阵信息查看前几行
> mi_df[1:4,1:4]
[,1] [,2] [,3] [,4]
miRNA_ID "hsa-let-7a-1" "hsa-let-7a-2" "hsa-let-7a-3" "hsa-let-7b"
read_count " 108784" " 109625" " 109658" " 148125"
miRNA_ID "hsa-let-7a-1" "hsa-let-7a-2" "hsa-let-7a-3" "hsa-let-7b"
read_count " 66736" " 66484" " 66714" " 64790"
这样的合并后得到的是有些冗余的结果，需要瘦身
> colnames(mi_df) <- mi_df[1,]
> #奇数列是多余的，只保留偶数列
> mi_df <- mi_df[seq(2,nrow(mi_df),2),]
> dim(mi_df)
[1] 116 1881
> mi_df[1:4,1:4]
hsa-let-7a-1 hsa-let-7a-2 hsa-let-7a-3 hsa-let-7b
read_count " 108784" " 109625" " 109658" " 148125"
read_count " 66736" " 66484" " 66714" " 64790"
read_count " 33740" " 33546" " 33772" " 34418"
read_count " 33770" " 33675" " 33476" " 19526"
> #转为数值型
> mi_df <- apply(mi_df, 2, as.numeric)
> mi_df[1:4,1:4]
hsa-let-7a-1 hsa-let-7a-2 hsa-let-7a-3 hsa-let-7b
[1,] 108784 109625 109658 148125
[2,] 66736 66484 66714 64790
[3,] 33740 33546 33772 34418
[4,] 33770 33675 33476 19526
最终的这个表达矩阵有点奇怪，就是没有样本信息，这个时候就需要对应到临床信息了。
继续学习：小洁老师讲可以通过TCGAbiolinks完成这个数据处理的过程，继续往后看。
