其他
TCGA2.GDC数据整理
大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~
就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~
这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！
花花写于2020.1.2 今天去看牙医了，痛苦啊。但是不耽误学习。
本系列是我的TCGA学习记录，跟着生信技能树B站课程学的，已获得授权（嗯，真的^_^）。课程链接：https://www.bilibili.com/video/av49363776
目录：TCGA-1.数据下载
1.xml文件探索
使用R包XML：
library(XML)
result <- xmlParse("./clinical/00e01b05-d939-49e6-b808-e1bab0d6773b/nationwidechildrens.org_clinical.TCGA-J2-8192.xml")
rootnode <- xmlRoot(result)
rootsize <- xmlSize(rootnode)
print(rootnode[2])
探索发现，xml共有两个节点，第二个节点中储存着病人的信息。
xmldataframe <- xmlToDataFrame(rootnode[2])
t(xmlToDataFrame(rootnode[2]))
[,1]
additional_studies ""
tumor_tissue_site "Lung"
histological_type "Lung Adenocarcinoma- Not Otherwise Specified (NOS)"
other_dx "Yes, History of Prior Malignancy"
gender "FEMALE"
vital_status "Alive"
2.写循环， 临床信息整理为数据框
xmls = dir("clinical/",pattern = "*.xml$",recursive = T)
td = function(x){
result <- xmlParse(file.path("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]
# additional_studies tumor_tissue_site histological_type
# [1,] "" "Lung" "Lung Adenocarcinoma- Not Otherwise Specified (NOS)"
# [2,] "" "Lung" "Lung Adenocarcinoma- Not Otherwise Specified (NOS)"
# [3,] "" "Lung" "Lung Adenocarcinoma- Not Otherwise Specified (NOS)"
3.探索miRNA信息
x = read.table(file = "mirna/00c46e8b-f303-4d20-bd6d-d650c36895f5/af29b644-b3a8-455b-9f07-b956d41f6ec4.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 75848 21277.9626 N
# 2 hsa-let-7a-2 75451 21166.5905 Y
# 3 hsa-let-7a-3 75769 21255.8004 N
# 4 hsa-let-7b 54045 15161.4741 N
# 5 hsa-let-7c 23828 6684.5704 Y
# 6 hsa-let-7d 1794 503.2785 N
第一列和第二列有用。探索一下会发现，每个病人的第一列miRNA_ID内容和顺序都是一样的。
4.写循环合成表达矩阵
mis = dir("mirna/",pattern = "*mirnas.quantification.txt$",recursive = T)
ex = function(x){
result <- read.table(file.path("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]
# [,1] [,2] [,3] [,4]
# miRNA_ID "hsa-let-7a-1" "hsa-let-7a-2" "hsa-let-7a-3" "hsa-let-7b"
# read_count " 43563" " 44020" " 43680" " 65476"
# miRNA_ID "hsa-let-7a-1" "hsa-let-7a-2" "hsa-let-7a-3" "hsa-let-7b"
# read_count " 75848" " 75451" " 75769" " 54045"
colnames(mi_df) <- mi_df[1,]
#奇数列是多余的，只保留偶数列
mi_df <- mi_df[seq(2,nrow(mi_df),2),]
dim(mi_df)
#[1] 567 1881
mi_df[1:4,1:4]
# hsa-let-7a-1 hsa-let-7a-2 hsa-let-7a-3 hsa-let-7b
# read_count " 43563" " 44020" " 43680" " 65476"
# read_count " 75848" " 75451" " 75769" " 54045"
# read_count " 56189" " 56527" " 56437" " 44829"
# read_count " 7956" " 7951" " 7834" " 15317"
#转为数值型
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,] 43563 44020 43680 65476
# [2,] 75848 75451 75769 54045
# [3,] 56189 56527 56437 44829
# [4,] 7956 7951 7834 15317
这样得到的表达矩阵有点问题，就是没有行名。这个问题应该可以得到解决，我的思路是用到列表的元素名称和TCGAid和下载的文件名的对应关系，在clinical信息中有，但有一点疑问，clinical和样本数量并不一致，给它增加行名可能丢失部分信息，所以先不解决，忍住~
总结
本文的结果是得到了clinical表格和miRNA表达矩阵。
下一篇我将分享用R包TCGAbiolinks实现这两步的办法。
插个小广告！
再给生信技能树打个call！