TCGA数据库：SNP数据的下载整理及其可视化
单核苷酸多态性主要是指在基因组水平上由单个核苷酸的变异所引起的DNA序列多态性。它是人类可遗传的变异中最常见的一种，占所有已知多态性的90%以上。SNP在人类基因组中广泛存在，平均每300个碱基对中就有1个，估计其总数可达300万个甚至更多。SNP是一种二态的标记，由单个碱基的转换或颠换所引起，也可由碱基的插入或缺失所致。SNP既可能在基因序列内，也可能在基因以外的非编码序列上。
SNP所表现的多态性只涉及到单个碱基的变异，这种变异可由单个碱基的转换(transition)或颠换(transversion)所引起，也可由碱基的插入或缺失所致。但通常所说的SNP并不包括后两种情况。
理论上讲，SNP既可能是二等位多态性，也可能是3个或4个等位多态性，但实际上，后两者非常少见，几乎可以忽略。因此，通常所说的SNP都是二等位多态性的。这种变异可能是转换(C←→T，在其互补链上则为G←→A)，也可能是颠换(C←→A，G←→T，C←→G，A←→T)。转换的发生率总是明显高于其它几种变异，具有转换型变异的SNP约占2/3，其它几种变异的发生几率相似。Wang等的研究也证明了这一点。转换的几率之所以高，可能是因为CpG二核苷酸上的胞嘧啶残基是人类基因组中最易发生突变的位点，其中大多数是甲基化的，可自发地脱去氨基而形成胸腺嘧啶。
在基因组DNA中，任何碱基均有可能发生变异，因此SNP既有可能在基因序列内，也有可能在基因以外的非编码序列上。总的来说，位于编码区内的SNP(coding SNP,cSNP)比较少，因为在外显子内，其变异率仅及周围序列的1/5。但它在遗传性疾病研究中却具有重要意义，因此cSNP的研究更受关注。
从对生物的遗传性状的影响上来看，cSNP又可分为2种：一种是同义cSNP(synonymous cSNP),即SNP所致的编码序列的改变并不影响其所翻译的蛋白质的氨基酸序列，突变碱基与未突变碱基的含义相同；另一种是非同义cSNP(non-synonymous cSNP),指碱基序列的改变可使以其为蓝本翻译的蛋白质序列发生改变，从而影响了蛋白质的功能。这种改变常是导致生物性状改变的直接原因。cSNP中约有一半为非同义cSNP。
先形成的SNP在人群中常有更高的频率，后形成的SNP所占的比率较低。各地各民族人群中特定SNP并非一定都存在，其所占比率也不尽相同，但大约有85%应是共通的。
本文介绍TCGA数据库中SNP的数据下载与整理
二.数据下载
这里下载方式采用网页筛选下载的方法，具体和之前的文章筛选条件一样。
不过，这里需要说明的是SNP的数据类型只能下载Masked Somatic Mutation，其他数据类型需要申请。所以我们这里下载的是Masked Somatic Mutation 的数据。
此外，在Workflow Type里面存在四个选项，MuSE Variant Aggregation and Masking，MuTect2 Variant Aggregation and Masking，SomaticSniper Variant Aggregation and Masking，VarScan2 Variant Aggregation and Masking，分别代表用四个软件计算出来的SNP的信息，四个软件分别采用四个不同的流程来处理数据计算得出的SNP信息，我们可以选择一种进行下载，比如我们选择VarScan2 Variant Aggregation and Masking，也可以下载所有的数据。
不过下载数据的时候，我们同时也要下载TSV格式的临床数据。因为整理SNP的数据需要结合临床数据一起整理。
关于临床数据的下载，参考文章：TCGA数据库：临床数据下载与整理。
下载后的数据解压会得到4个文件夹，也就是上面提到的4个软件分析的数据。我们选择其中一个(VarScan2 Variant Aggregation and Masking)进行讲解。
三.数据的整理
1.设置工作路径 我们下载的数据是个压缩文件，解压后是maf的格式文件。所以注意文件的输入与输出路径。
rm(list = ls())
options(stringsAsFactors = F)
setwd("F:\\TCGAcourse\\01-dataDownload\\SNP")2.读入临床数据
选择其中一部分进行演示，更多处理参考：TCGA数据库：临床数据下载与整理，也要注意临床数据的路径。
# 处理临床数据
clin = read.table(".\\COAD\\clinical.cart\\clinical.tsv",header = T,sep = "\t")
clin1 <- clin[,c("case_id","submitter_id","vital_status","days_to_last_follow_up",
"ajcc_pathologic_m","ajcc_pathologic_n",
"ajcc_pathologic_stage","ajcc_pathologic_t")]
colnames(clin1) = c("case_id","Tumor_Sample_Barcode","vital_status","days_to_last_follow_up",
"ajcc_pathologic_m","ajcc_pathologic_n",
"ajcc_pathologic_stage","ajcc_pathologic_t")
clin1 <- clin1[clin1$days_to_last_follow_up!= "--",]
clin1$days_to_last_follow_up <- as.numeric(clin1$days_to_last_follow_up)
clin1 <- clin1[clin1$vital_status!= "--",]
clin1$vital_status <- ifelse(clin1$vital_status== "Alive",0,1)3.读入数据
读入数据我们需要maftools包中的read.maf函数。clinicalData就是上一步读入的临床数据
library(maftools)
laml = read.maf(maf = ".\\COAD\\TCGA.COAD.varscan.somatic.maf",
clinicalData = clin1,
isTCGA = TRUE)读入的数据是一个列表类型的数据。
我们可以提取SNP中病人对应的临床数据：
clindata = laml@clinical.data ###提取SNP中对应病人的临床数据可以通过plotmafSummary函数查看整体统计摘要。
plotmafSummary(laml, rmOutlier = TRUE, dashboard = TRUE,
titvRaw = TRUE, addStat = NULL, showBarcodes = FALSE, fs = 0.8,
textSize = 0.8, color = NULL, titleSize = c(0.8, 0.6),
titvColor = NULL, top = 10)
也可以提取特定基因的突变信息：
# 提取肿瘤样本中特定基因的SNP统计数据
tp53 = genesToBarcodes(laml, genes = "TP53")[[1]]
下面这些函数，运行一下就知道是什么意思啦。
#使用等位基因频率或按突变状态创建基因型矩阵。
KRAS = genotypeMatrix(maf = laml, genes = "KRAS")
# 从MAF对象中提取注释,其实就是Barcode
Clin = getClinicalData(laml)
#其实就是lama对象中data的列名
colNM = getFields(laml)
##
mafCompare(laml,laml)
###总结基因和样本，不考虑改变的类型。
mafs <- mafSummary(laml)
#基因突变数据摘要
lamaGene = getGeneSummary(laml)
#样本突变数据摘要
lamaSample = getSampleSummary(laml)可能重要的是我们需要提取数据
snpdata <- laml@data
关于列名，看名称就能看的懂，Variant_Type这一列就突变类型。
4.数据整理
我们只提取下面4列
tidydata <- snpdata[,c("Hugo_Symbol","Variant_Type",
"Variant_Classification","Tumor_Sample_Barcode")]删除无义突变，因为这个没有意义。
#删除无义突变
tidydata <- tidydata[which(tidydata$Variant_Type != "Nonsense_Mutation"),
c("Hugo_Symbol","Variant_Type","Tumor_Sample_Barcode")]然后我们利用reshape2包对数据进行重塑。0代表野生型，其他数字代表有突变，1表示该基因在该病人中有1个突变，2代表有2个突变。
library(reshape2)
rsdata = dcast(tidydata, Hugo_Symbol~Tumor_Sample_Barcode)
rownames(rsdata) <- rsdata$Hugo_Symbol
rsdata <- rsdata[,-1] #0代表野生型，其他数字代表有突变
得到这样的数据后，我们就可以进行后续的分析，比如结合临床数据的分析。也可以结合RNA-Seq进行分析等。后续再慢慢介绍吧。另外，TCGA数据库33个Project的RNA-Seq转录组数据为你整理打包好了，需要的可以下载，处理方式参考文章：TCGA数据库：RNA-Seq数据的下载与处理，差异分析，参考文章：一文就会TCGA数据库基因表达差异分析。
本文分享自 MedBioInfoCloud 微信公众号，前往查看
如有侵权，请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体分享计划 ，欢迎热爱写作的你一起参与！