24

一个标准的TCGA大文章应该做哪些数据?

很多人总是问我如何挖掘TCGA的数据,发文章!
可是他却连TCGA的数据是怎么来的都不知道,TCGA发了几十篇CNS大文章(自己测序的)了,每篇文章都有几百个左右的癌症样本的6种数据,这几年凑成了一万多个样本,都放在GDC里面可以任意下载。同时也出来了十几篇TCGA的数据挖掘大文章(主要包括亚型,driver mutation,假基因等新型研究领域)
那么一篇标准的一个标准的TCGA大文章应该自己测哪些数据?
其实稍微仔细浏览几篇文章就明白了,套路也是存在的,https://tcga-data.nci.nih.gov/docs/publications/
我们就以2013年发表在新英格兰杂志上面的Genomic and Epigenomic Landscapes of Adult De Novo Acute Myeloid Leukemia 为例子吧!

Continue reading

十二 28

TCGA表达数据的多项应用之4–求指定基因在指定癌症里面的表达量相关性矩阵,与所有的基因比较。

这个不出图,会给出TCGA里面涉及到的所有基因跟你指定的基因的表达量相关系数和P值,分别你一次性的看清楚你感兴趣的基因跟体内其它基因在该癌症种类的相关性,当然,相关非因果,请谨慎应用! Continue reading

十二 25

TCGA表达数据的多项应用之1–下载数据并且导入mysql

这个TCGA表达数据的多项应用系列帖子是应群里朋友的要求来写的,你们也可以继续提需求,我会接着写下去,其实从TCGA数据库里面下载到了数据之后,后面的所有分析都跟TCGA没有半毛钱关系了,大家要有这个想法,别三两句就问TCGA数据怎么分析,http://www.bio-info-trainee.com/?s=TCGA&submit=Search 本系列最后会形成一个shiny版本的交互式表达数据查询,处理,绘图,统计的网页APP。
我这里偷懒一下了,直接下载GEO里面的TCGA的表达数据,而不是去TCGA的官网里面下载:
它处理了目前(大概是2015年6月)TCGA收集的所有癌症样本的mRNA表达数据,并且统一处理成了count和RPKM两种表达量形式。 GEO地址:http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62944

Continue reading

19

用GISTIC多个segment文件来找SCNA变异

这个软件在TCGA计划里面被频繁使用者,用这个软件的目的很简单,就是你研究了很多癌症样本,通过芯片得到了每个样本的拷贝数变化信息,芯片结果一般是segment结果,可以解释为CNV区域,需要用GISTIC把样本综合起来分析,寻找somatic的CNV,并且注释基因信息。

有两个难点,一是在linux下面安装matlab工作环境,二是如何制作输入文件。

Continue reading

18

2016-TCGA数据挖掘系列文章之癌症男女有别

这是TCGA数据挖掘系列文章之一,是安德森癌症研究中心的Han Liang主导的,纯粹的生物信息学数据分析文章。
文章题目是:comprehensive characterization of molecular differences in cancer between male and female patients.
研究意义:癌症病人的性别对肿瘤发生,扩散的意义不言而喻。不仅仅是因为很多癌症本来就是有性别特异性,比如卵巢癌之于女性、前列腺癌之于男性。即使对于其它并非性别特异性的癌症种类,男女病人在肿瘤发生,扩散,以及治疗阶段的反应也大不一样。但是以前对这样分子机理研究的很有限,一般集中在某些性别相关的分子pattern,比如非小细胞肺癌女性患者的EGFR突变,但那些研究要么就局限于单一的基因,要么局限于单一的数据类型,或者研究单一的癌症。严重缺乏一个全面的,系统的分析癌症患者的性别差异。而且TCGA数据库的出现让这一个研究变成了可能,这也就是本文章的出现的原因。
数据挖掘的对象:
如表所示,涉及到13种癌症,TCGA的六种数据()都用上了,因为是2016年,所以数据量也比较全面了。

Continue reading

16

TCGA数据挖掘系列文章之-pseudogene假基因探究

这是TCGA数据挖掘系列文章之一,是安德森癌症研究中心的Han Liang主导的,纯粹的生物信息学数据分析文章。
TCGA数据库的数据量现在已经非常可观了,一万多的肿瘤样本数据,关于假基因的这篇文章是2014年发的,所以他们只研究了2,808个样本数据,也只涉及到7个癌症种类。
假基因是原来的能翻译成蛋白的基因经过各种突变导致丧失功能的基因。
比如
PTEN-->PTENP1
KRAS-->KRASP1
NANOG-->NANOGP1
很好理解,一般来说看到结尾是P1,等字眼的都是假基因,现在共有一万多假基因,我一般以http://www.genenames.org/cgi-bin/statistics (人类基因命名委员会)为标准参考。
文章主要做了6件事情。
一是重新定义及规范了假基因该研究什么就是把Yale Pseudogene database的假基因资源和GENCODE Pseudogene Resource的假基因资源结合起来,然后定义了一些过滤手段,具体流程如下。
1
二是下载了TCGA的那2,808个样本的RNA-seq的level2数据,也就是bam文件,重新提取关于假基因的表达数据。如果只是自己下载表达数据的话,关于假基因的定量并不准确,而且只有五百多个假基因。
当然,一般人没有条件下载RNA-seq的level2数据,所以想学习这个流程的话,直接下载表达矩阵吧。
Cancer type Number of nontumour samples Number of tumour samples Sequencing strategy Number of mappable reads Number of detectable pseudogenes
Breast invasive carcinoma 105 837 Paired-end 161 M 747
Kidney renal clear cell carcinoma 67 448 Paired-end 166 M 712
Lung squamous cell carcinoma 17 220 Paired-end 171 M 813
Ovarian serous cystadenocarcinoma 0 412 Paired-end 170 M 670
Glioblastoma multiforme 0 154 Paired-end 106 M 875
Colorectal carcinoma 0 228 Single-end 22 M 168
Uterine corpus endometrioid carcinoma 4 316 Single-end 26 M 181
第三件事是把假基因与其配对的野生型基因的表达数据做了相关性分析,一般来说,它们的相关性由下面三个原因决定。
(i) the sequence similarity between the pseudogene/gene pair;
(ii) the molecular mechanisms through which the pseudogene functions;
(iii) the detection sensitivity given the setting of RNA-seq experiments.
结论是不怎么相关,暗示着假基因虽然不编码蛋白产物,但仍然行使着某种功能。
第四件事是如果RNA-seq有正常对照的, 就做一样normal和tumor的差异分析,当然现在已经是都有了,在GSE62944可以下载所有的表达数据,专门提取假基因的表达数据做差异分析就好了。
但是差异分析的结果是, 没有什么现实意义。所以作者认为normal和tumor这样比较是不科学的,因为tumor本来就不应该按照组织来分类,而是应该按照TCGA的6种数据来分类()
In recent years, various ‘omic’ data, such as mRNA expression, microRNA expression, DNA methylation, somatic copy number alteration and protein expression, have been widely used to classify tumour samples into different molecular subtypes13, 14, 15, 16, 17, 18, 19.
2
第五件事就是把假基因表达数据的分类来跟其它几种分类形式作比较。
那些分类来源于以前的TCGA大文章:
48 in UCEC (endometrioid vs serous)23,
138 in LUSC (basal, classical, primitive and secretory)16,
71 in GBM (classical, mesenchymal, neural and proneural)24 and
547 in BRCA (PAM50 subtypes: luminal A, luminal B, basal-like, Her2-enriched and normal-like)2
文章就是:13141516171819.
3
最后就是做一些生存分析,讲一些好听的故事,比如说这样分类有利于精准医疗。

看起来还不错,值得大家学习一下,数据也都可以下载, 文章中提供了syn编号。

06

所有TCGA的maf格式somatic突变数据均可下载

如果你研究癌症,那么TCGA计划的如此丰富的公共数据你肯定不能错过,一般人只能获取到level3的数据,当然,其实一般人也没办法使用level1和level2的数据,毕竟近万个癌症样本的原始测序数据,还是很恐怖的,而且我们拿到原始数据,再重新跑pipeline,其实并不一定比人家TCGA本身分析的要好,所以我们直接拿到分析结果,就足够啦!

Continue reading

06

突变频谱探究mutation siganures

这也是对TCGA数据的深度挖掘,从而提出的一个统计学概念。文章研究了30种癌症,发现21种不同的mutation signature。如果理解了,就会发现这个其实蛮简单的,他们并不重新测序,只是拿已经有了的TCGA数据进行分析,而且居然是发表在nature上面!

研究了4,938,362 mutations from 7,042 cancers样本,突变频谱的概念只是针对于somatic 的mutation。一般是对癌症病人的肿瘤组织和癌旁组织配对测序,过滤得到的somatic mutation,一般一个样本也就几百个somatic 的mutation。

paper链接是:http://www.nature.com/nature/journal/v500/n7463/full/nature12477.html

从2013年提出到现在,已经有30种mutation siganures,在cosmic数据库有详细记录,更新见:http://cancer.sanger.ac.uk/cosmic/signatures
它的概念就是:根据突变上下文分成96类,然后每类突变的频率不一样画一个条形图,可视化展现。
mutation signature

Each signature is displayed according to the 96 substitution classification defined by the substitution class and sequence context immediately 3′ and 5′ to the mutated base. The probability bars for the six types of substitutions are displayed in different colours.
仔细看paper,还是蛮好理解的,自己写一个脚本就可以做这个分析了,前提是下载各个癌症的somatic mutation文件,一般是maf格式的,很多途径下载。
In principle, all classes of mutation (such as substitutions, indels, rearrangements) and any accessory mutation characteristic, for example, the sequence context of the mutation or the transcriptional strand on which it occurs, can be incorporated into the set of features by which a mutational signature is defined. In the first instance, we extracted mutational signatures using base substitutions and additionally included information on the sequence context of each mutation. Because there are six classes of base substitution—C>A, C>G, C>T, T>A, T>C, T>G (all substitutions are referred to by the pyrimidine of the mutated Watson–Crick base pair)—and as we incorporated information on the bases immediately 5′ and 3′ to each mutated base, there are 96 possible mutations in this classification. This 96 substitution classification is particularly useful for distinguishing mutational signatures that cause the same substitutions but in different sequence contexts.

很多癌症都发现了不止一种mutation signature,甚至高达6种,说明癌症之间差异还是蛮大的!
In most cancer classes at least two mutational signatures were observed, with a maximum of six in cancers of the liver, uterus and stomach. Although these differences may, in part, be attributable to differences in the power to extract signatures, it seems likely that some cancers have a more complex repertoire of mutational processes than others.
Most individual cancer genomes exhibit more than one mutational signature and many different combinations of signatures were observed
但是,我最后也没能绝对的界限是什么,因为总不能用肉眼来看每个突变频谱不一样吧?
The set of signatures will be updated in the future. This will include incorporating additional mutation types (e.g., indels, structural rearrangements, and localized hypermutation such as kataegis) and cancer samples. With more cancer genome sequences and the additional statistical power this will bring, new signatures may be found, the profiles of current signatures may be further refined, signatures may split into component signatures and signatures may be found in cancer types in which they are currently not detected.
分类会持续不断更新,随着更多的cancer type和样本加入,新的signature会被发现,现有的signature也可能会被重新定义,或者被分割成多个小的signature
28

使用可视化工具MutationMapper来看看基因上面突变的分布

how to generate lollipop diagrams ?
这个工具是网页版的,不需要下载,只需上传两列数据即可,就是基因上面的第几个氨基酸是如何突变的这样的信息
网页就可以把这些信息画在基因上面,而且注释到domain和cosmic数据库,挺好用的!
Hugo_Symbol HUGO symbol for the gene TP53
Protein_Change Amino acid change V600E
结果也很容易理解!
我输入了4个基因,网页就会对每个基因画一个图,其实ABCA1这个基因有6个突变,就都画在该基因上面了,基因突变分布在基因的哪个domain也看得一清二楚!
1
我的输入文件是这样的!
2
非常好用:
  • Support mutation data with annotated protein effects
  • Mutation diagram/lollipop view
  • Mutation table view
  • 3D structure view if available
而且pfam数据库本身也提供了这个功能:

Pfam provides an online tool to not only generate the domain information in JSON format, but to draw the lollipop diagram using javascript as well. They have more information here: http://pfam.xfam.org/help#tabview=tab9

IMHO, not as pretty as cBioPortal's but it gets you close to a solution.

EDIT / SHAMELESS PLUG: After seeing the data available and how easy it'd be, I made my own quick tool to fetch the data and draw the diagram for me in a style similar to cBioPortal - feel free to fork it and add features: https://github.com/pbnjay/lollipops

Example output (w/ labels per the comments)

3

甚至还有高手用D3.JS写了一个能实现同样需求的模块

We found ourselves in the same need, we wanted such a plot (JavaScript). Thus, I add our solution, Mutations Needle Plot. The library creates an SVG image (with D3), which then may be downloaded.

You will npm in order to be able to install & run the library.

Examples may be found in the snippets folder or also the index.html - The one displayed here below

类似的高手很多:
My colleague @SolenaLS recently asked me to write something like this: (uniprot+SVG+javascript : )http://lindenb.github.io/pages/uniprot/paintsvg.html
而且ensembl数据库也提供类似的功能
Ensemble also gives similar plot.

 

22

用TCGA数据做cox生存分析的风险因子(比例风险模型)

再次强调一下,R里面实现生存分析非常简单!
用my.surv <- surv(OS_MONTHS,OS_STATUS=='LIVING')构建生存曲线。
用kmfit2 <- survfit(my.surv~TUMOR_STAGE_2009)来做某一个因子的KM生存曲线。用 survdiff(my.surv~type, data=dat)来看看这个因子的不同水平是否有显著差异,其中默认用是的logrank test 方法。
用coxph(Surv(time, status) ~ ph.ecog + tt(age), data=lung) 来检测自己感兴趣的因子是否受其它因子(age,gender等等)的影响。

Continue reading

14

用RNA-SeQC得到表达矩阵RPKM值

这个软件不仅仅能做QC,而且可以统计各个基因的RPKM值!尤其是TCGA计划里面的都是用它算的
一、程序安装
直接在官网下载java版本软件即可使用:http://www.broadinstitute.org/cancer/cga/tools/rnaseqc/RNA-SeQC_v1.1.8.jar
但是需要下载很多注释数据
clipboard

二、输入数据

clipboard

箭头所指的文件,一个都不少,只有那个rRNA.tar我没有用, 因为这个软件有两种使用方式,我用的是第一种
三、软件使用
软件的官网给力例子,很容易学习:
RNA-SeQC can be run with or without a BWA-based rRNA level estimation mode. To run without (less accurate, but faster) use the command:
java -jar RNASeQC.jar -n 1000 -s "TestId|ThousandReads.bam|TestDesc" -t gencode.v7.annotation_goodContig.gtf -r Homo_sapiens_assembly19.fasta -o ./testReport/ -strat gc -gc gencode.v7.gc.txt 
我用的就是这个例子,这个例子需要的所有文件里面,染色体都是没有chr的,这个非常重要!!!
代码如下:
 java -jar RNA-SeQC_v1.1.8.jar  \
-n 1000 \
-s "TestId|ThousandReads.bam|TestDesc" \
-t gencode.v7.annotation_goodContig.gtf \
-r ~/ref-database/human_g1k_v37/human_g1k_v37.fasta  \
-o ./testReport/ \
-strat gc \
-gc gencode.v7.gc.txt \
To run the more accurate but slower, BWA-based method :
java -jar RNASeQC.jar -n 1000 -s "TestId|ThousandReads.bam|TestDesc" -t gencode.v7.annotation_goodContig.gtf -r Homo_sapiens_assembly19.fasta -o ./testReport/ -strat gc -gc gencode.v7.gc.txt -BWArRNA human_all_rRNA.fasta
Note: this assumes BWA is in your PATH. If this is not the case, use the -bwa flag to specify the path to BWA
四、结果解读
运行要点时间,就那个一千条reads的测试数据都搞了10分钟!
出来一大堆突变,具体解释,官网上面很详细,不过,比较重要的当然是RPKM值咯,还有QC的信息
clipboard
TCGA数据里面都会提供由RNA-SeQC软件处理得到的表达矩阵!
Expression
  • RPKM data are used as produced by RNA-SeQC.
  • Filter on >=10 individuals with >0.1 RPKM and raw read counts greater than 6.
  • Quantile normalization was performed within each tissue to bring the expression profile of each sample onto the same scale.
  • To protect from outliers, inverse quantile normalization was performed for each gene, mapping each set of expression values to a standard normal.
软件的主页是:
 
 
 
 
11

对CCLE数据库可以做的分析

收集了那么多的癌症细胞系的表达数据,拷贝数变异数据,突变数据,总不能放着让它发霉吧!
这些数据可以利用的地方非常多,但是在谷歌里面搜索引用了它的文章却不多,我挑了其中几个,解读了一下别人是如何利用这个数据的,当然,主要是用那个mRNA的表达数据咯!
这篇文献对CCLE的数据进行了八个步骤的处理,一个合格的生物信息学分析着完全可以重写这个过程
step1:Affymetrix U133 Plus2 DNA microarray gene expressions of 27 gastric cancer cell lines (Kato-III, IM95, SNU-620, SNU-16, OCUM-1, NUGC-4, 2313287, HUG1N, MKN45, NCIN87, KE39, AGS, SNU-5, SNU-216, NUGC-3, NUGC-2, MKN74, MKN7, RERFGC1B, GCIY, KE97, Fu97, SH10TC, MKN1, SNU-1, Hs746 T, HGC27) were downloaded from Cancer Cell Line Encyclopedia (CCLE) [16] in March 2013.
step2: Robust Multi-array Average (RMA) normalization was performed. Principal component analysis plot show no obvious batch effect.
step3: The normalized data is then collapsed by taking the probe sets with highest gene expression.
前三步是为了得到27个胃癌相关细胞系的mRNA表达矩阵,方法是下载cel文件,用RMA归一化,对多探针基因去最大表达量探针!

step4:Unsupervised hierarchical clustering (1-Spearman distance, average linkage) was performed on the cell lines using the aCGH data.

Putative driver genes of which copy number aberrations correlated to mRNA gene expression were identified to determine subtypes or clusters that are driven by different mechanisms. This was done using Mann Whitney U-test with p<0.05, and Spearman Correlation Coefficient test with Rho >0.6.

step5:We then performed consensus clustering[17] on the gene expression data of the 27 gastric cancer cell lines from CCLE using these putative driver genes. We selected k = 2 as it gives sufficiently stable similarity matrix.

step6: In order to assign new samples to this integrative cluster, significance analysis of microarray (SAM) [18]with threshold q<2.0 was used to generate subtype signature based on the mRNA expression data of the 1762 genes from the 27 gastric cancer cell lines in CCLE.

先用甲基化数据来聚类,得到putative driver genes,然后再用这些基因的表达数据来再次聚类,分成两类,然后对这两类进行SAM找差异基因

step7:ssGSEA (single sample GSEA)was used to estimate pathway activities of the gastric cancer cell line in the Molecular Signature Database v3.1 (Msigdb v3.1) [19][20]. The pathway activities are represented in enrichment scores which were rank normalized to [0.0, 1.0]. 
step8:SAM analysis was performed with threshold q<0.2, and fold change >2.0 (for up-regulated pathways), or <0.5 (for down-regulated pathways) to obtain subtype-specific pathways from the 27 gastric cell lines in CCLE.
这里既用来gene set的富集分析,又用来超几何分布的富集分析,结果去看看这篇文章就知道了!
 
这篇文章只用了CCLE的一个地方,就是看看不同cancer type里面的某个基因表达boxplot
这个图的数据用GEOquery可以得到,样本的分类信息也用GEOquery可以得到,这样就可以做下面这个图了,非常简单
Further, the Cancer Cell Line Encyclopedia (CCLE) database demonstrated that of 1062 cell lines representing 37 distinct cancer types, glioma cell lines express the highest levels of STK17A
1

结论就是:STK17A is highly expressed in glioma cell lines compared to other cancer types. Data was obtained through the Cancer Cell Line Encyclopedia (CCLE).

第三篇文献:http://www.nature.com/ncomms/2013/130709/ncomms3126/fig_tab/ncomms3126_F4.html

这篇文献更简单了,直接对这个表达矩阵进行聚类:
 
The 5,000 most variable genes were used for unsupervised clustering of cell lines by mRNA expression data. Cell lines are colour-coded (vertical bars) according to the reported tissue of origin (a PDF version that can be enlarged at high resolution is in Supplementary InformationSupplementary Fig. S4); horizontal labels at bottom indicate the dominating tissue types within the respective branches of the dendrogram. Most ovarian cancer cell lines (magenta) cluster together, interspersed with endometrial cell lines. However, some ovarian cancer cell lines cluster with other tissue types (*). Top right panels: neighbourhoods (1) of the top cell lines in our analysis, (2) of cell line IGROV1, and (3) of cell line A2780. For the ovarian cancer cell lines in these enlarged areas, the histological subtype as assigned in the original publication is indicated by coloured letters.
就直接拿整个表达矩阵即可,然后挑选变异最大的5000个基因来进行聚类,就可以得到类似的图

 

11

CCLE数据库几个知识点

Here we describe the Cancer Cell Line Encyclopedia (CCLE): a compilation of gene expression, chromosomal copy number, and massively parallel sequencing data from 947 human cancer cell lines. 
收集了三种数据:
The mutational status of >1,600 genes was determined by targeted massively parallel sequencing, followed by removal of variants likely to be germline events . 
Moreover, 392 recurrent mutations affecting 33 known cancer genes were assessed by mass spectrometric genotyping13 . 
DNA copy number was measured using high-density single nucleotide polymorphism arrays (Affymetrix SNP 6.0; Supplementary Methods). 
Finally, mRNA expression levels were obtained for each of the lines using Affymetrix U133 plus 2.0 arrays. 
These data were also used to confirm cell line identities .
一般用得最多的就是表达数据,因为表达数据最简单,大多数生物信息学分析着只会用这个数据!
而它的突变数据又不是通常意义的高通量测序得到的,snp6芯片数据很多人听都没听过
文章的附件有对cell lines的具体描述。
different_kinds_of_cancer_in_CCLE
CCLE的数据在broad institute里面可以下载,也放在GEO数据库里面,我比较喜欢GEO里面的数据
This SuperSeries is composed of the following SubSeries:
GSE36133 Expression data from the Cancer Cell Line Encyclopedia (CCLE)
GSE36138 SNP array data from the Cancer Cell Line Encyclopedia (CCLE)
GSE36133这个study的metadata里面有对每个cellline来源的cancer进行描述!
有人喜欢把这个metadata叫做是clinical data。
library(GEOquery)
ccleFromGEO <- getGEO("GSE36133")
annotBlock1 <- pData(phenoData(ccleFromGEO[[1]]))
>dim(annotBlock1)
[1] 917  38
exprSet=exprs(ccleFromGEO[[1]])
> dim(exprSet)
[1] 18926   917
##它的表达数据矩阵,包含了18926个基因,列名是917个细胞系的名字,行是基因的entrez ID
keyColumns <- c("title", "source_name_ch1", "characteristics_ch1", "characteristics_ch1.1", 
    "characteristics_ch1.2")
options(stringsAsFactors = F)
allAnnot=annotBlock1[,keyColumns]
##这几列信息是比较重要的metadata,里面详细记录了细胞系的收集公司单位,tissue,癌症分类等信息
Cell line (1035个细胞系简介)Gene Sets
1035 sets of genes with high or low expression in each cell line relative to other cell lines from the CCLE Cell Line Gene Expression Profiles dataset.
一些关于CCLE数据库的文章:
http://onlinelibrary.wiley.com/doi/10.1002/cncy.21471/pdf 介绍了几个类似的数据库资源
Anticancer drug sensitivity analysis: An integrated approach applied to Erlotinib sensitivity prediction in the CCLE database
08

TCGA数据里面的生存分析例子

我们知道了生存分析,就是随着时间的流逝,死亡率是如何增加的,一般是用KM法来估计生存函数,然后画个图即可!而根据某些因子把样本分组,可以看到他们死亡率的变化趋势显著的不同,这就说明了我们的这个因子是非常有效的分类方式,这个因子可以是一个biomarker,也可以某些其它指标!
甚至,我们还可以用cox模型来分析这个因子是如何影响生存函数的,那个稍后再讲
这里,我们就简单讲一个例子,是TCGA里面卵巢癌的数据,根据甲基化数据分成了4个组,那么我们就下载这四个组样本的临床数据,
看看这样分组后,他们的死亡率变化趋势是不是有显著区别!
数据是用cgdsr下载的:
这个包的说明见:http://www.bio-info-trainee.com/?p=1257 

library(cgdsr)
test(mycgds)
all_TCGA_studies <- getCancerStudies(mycgds)
all_tables <- getCaseLists(mycgds, 'ov_tcga_pub')
all_dataset<- getGeneticProfiles(mycgds, 'ov_tcga_pub')
#BRCA1 <- getProfileData(mycgds, my_gene, my_dataset, my_table)
ov_tcga_pub_meth1<- getClinicalData(mycgds, all_tables[8,1])
ov_tcga_pub_meth2<- getClinicalData(mycgds, all_tables[9,1])
ov_tcga_pub_meth3<- getClinicalData(mycgds, all_tables[10,1])
ov_tcga_pub_meth4<- getClinicalData(mycgds, all_tables[11,1])
下载之后的数据如下:
1
根据甲基化数据,把癌症病人分成了4组,我们的临床数据记录了13项,但是我们只需要用到OS_MONTHS和OS_STATUS就可以来估计KM生存函数,画出生存曲线啦!
无病生存期(Disease-free  survival,DFS)的定义是指从随机化开始至疾病复发或由于疾病进展导致患者死亡的时间。该指标也常作 为抗肿瘤药物III期临床试验的主要终点。某些情况下,DFS与OS相比,作为终点比较难以记录,因为它要求认真随访,及时发现疾病复发,而且肿瘤患者的 死亡原因也很难确定。肿瘤患者常有合并症(如,心血管病),这些合并症可能会干扰对DFS的判断。并且,肿瘤患者常死于医院外,不能常规进行尸检。
总生存期(Overall survival,OS)的定义是指从随机化开始至因任何原因引起死亡的时间。该指标常常被认为是肿瘤临床试验中最佳的疗效终点。如果在生存期上有小幅度的提高,可以认为是有意义的临床受益证据。作为一个终点,生存期应每天进行评价,可通过在住院就诊时,通过与患者直接接触或者通过电话与患者交谈,这些相对比较容易 记录。确认死亡的日期通常几乎没有困难,并且死亡的时间有其独立的因果关系。当记录至死亡之前的失访患者,通常截止到最后一次有记录的、与患者接触的时间。
library(survival)
attach(ov_tcga_pub_meth1)
## 估计KM生存曲线
my.surv <- Surv(OS_MONTHS,OS_STATUS=='LIVING')
kmfit1 <- survfit(my.surv~1)
summary(kmfit1)
plot(kmfit1)
2
##我们很容易看到,随着感染癌症的时间延长,病人的死亡率到了一定程度就不增加了,有些病人熬过了这个癌症,或者说,到我们拿到数据为止,他们还没有死亡!
## 随便取一根因子来分组TUMOR_STAGE_2009估计KM生存曲线
kmfit2 <- survfit(my.surv~TUMOR_STAGE_2009)
summary(kmfit2)
plot(kmfit2,col = rainbow(length(unique(ov_tcga_pub_meth1[,13]))))
detach(ov_tcga_pub_meth1)
3
##可以看到,我们根据病人的TUMOR_STAGE_2009把他们分成了这些组,不同组的生存曲线不一样,但是我们肉眼无法看出它们这些组直接的生存率是否有显著差异!我们需要做统计检验!
我们就不对这个进行检验了,我们还是用下面的代码来对甲基化数据的分组来做检验,看看我们的分组是否有效!
ov_tcga_pub_meth1$sample<- rownames(ov_tcga_pub_meth1)
ov_tcga_pub_meth2$sample<- rownames(ov_tcga_pub_meth2)
ov_tcga_pub_meth3$sample<- rownames(ov_tcga_pub_meth3)
ov_tcga_pub_meth4$sample<- rownames(ov_tcga_pub_meth4)
ov_tcga_pub_meth1$type<- 'meth1'
ov_tcga_pub_meth2$type<- 'meth2'
ov_tcga_pub_meth3$type<- 'meth3'
ov_tcga_pub_meth4$type<- 'meth4'
dat=rbind(ov_tcga_pub_meth1,ov_tcga_pub_meth2,ov_tcga_pub_meth3,ov_tcga_pub_meth4)
attach(dat)
## 根据分组type估计KM生存曲线
my.surv <- Surv(OS_MONTHS,OS_STATUS=='LIVING')
kmfit3 <- survfit(my.surv~type)
summary(kmfit3)
plot(kmfit3,col = rainbow(4))
4
##由图中可以看到甲基化数据分成的4个组,生存率差异还是蛮大的
##用图形方法检验PH假设
plot(kmfit3,fun='cloglog',col = rainbow(4))
# 检验显著性
survdiff(my.surv~type, data=dat)
# 用strata来控制协变量的影响
survdiff(my.surv~type+strata(TUMOR_STAGE_2009),data=dat)
5
##然后我们用这个代码来检验,是否显著,结果发现,p值并没有小于0.05,只能说是可以这样分组,但是分组的效果没有我们想象中的那么好!
names(dat)
detach(dat)
08

生存分析简介

一般我们谈生存分析,就是说的KM方法估计生存函数,并且画出生存曲线,然后还可以根据分组检验一下它们的生存曲线是否有显著的差异!
在R里面,非常的方便,有个包survival很容易就可以做生存分析了!
只需要记住三个函数即可:
Surv:用于创建生存数据对象
survfit:创建KM生存曲线或是Cox调整生存曲线
survdiff:用于不同组的统计检验
我们用下面的数据做例子
clipboard
对于上面的数据,我们用下面的代码做生存分析!
library(survival)
my.surv <- Surv(OS_MONTHS,OS_STATUS=='LIVING')  ##这个生存对象是看看病人的总生存期死亡状态的关系
###这个Surv函数第一个参数必须是数值型的时间,第二个参数是逻辑向量,1,0表示死亡与否
kmfit1 <- survfit(my.surv~1) ### 直接对生存对象拟合生存函数
summary(kmfit1)
plot(kmfit1)   ###画出生存曲线

survdiff(my.surv~type, data=dat)     ### 根据生存对象再加上一个分组因子来拟合生存函数,并且比较不同因子分组的生存效果

 

05

寻找somatic突变的软件大合集

         其实somatic突变很容易理解,你测同一个人的正常组织和癌症组织,然后比较对这两个样本数据call出来的snp位点
         只存在癌症组织数据里面的snp位点就是somatic突变,在两个样本都存在的snp位点就是germline的突变,不过一般大家研究的都是somatic突变。
          当然,理论上是很简单,但是那么多统计学家要吃饭呀,肯定得把这件事搞复杂才行,所以就有了非常多的somatic突变 calling的软件,开个玩笑哈,主要是因为我们的测序并不是对单个细胞测序,我们通常意义取到的正常组织和癌症组织都不是纯的,所以会有很多关于这一点的讨论。
        正好我看到了一篇帖子,收集了大部分比较出名的做somatic mutation calling的软件,当然,我只用过mutect和varscan。

来自于:https://www.biostars.org/p/19104/

Here are a few more, a summary of the other answers, and updated links:

For a much more general discussion of variant calling (not necessarily somatic or limited to SNVs/InDels) check out this thread: What Methods Do You Use For In/Del/Snp Calling?

Some papers describing comparisons of these callers:

The ICGC-TCGA DREAM Mutation Calling challenge has a component on somatic SNV calling.

This paper used validation data to compare popular somatic SNV callers:

Detecting somatic point mutations in cancer genome sequencing data: a comparison of mutation callers

You'll need to update the link to MuTect. Broad Institute has begun to put portable versions of their tools on Github, like thelatest release of MuTectThe Genome Institute at WashU has been using Github for a while, but portable versions of their tools can be found here and here.

其实somatic的calling远比我们想象的要复杂:

To rehash/expand on what Dan said, if you're sequencing normal tissue, you generally expect to see single-nucleotide variant sites fall into one of three bins: 0%, 50%, or 100%, depending on whether they're heterozygous or homozygous.

With tumors, you have to deal with a whole host of other factors:

  1. Normal admixture in the tumor sample: lowers variant allele fraction (VAF)
  2. Tumor admixture in the normal - this occurs when adjacent normals are used, or in hematological cancers, when there is some blood in the skin normal sample
  3. Subclonal variants, which may occur in any fraction of the cells, meaning that your het-site VAF might be anywhere from 50% down to sub-1%, depending on the tumor's clonal architecture and the sensitivity of your method
  4. Copy number variants, cn-neutral loss of heterozygosity, or ploidy changes, all of which again shift the expected distribution of variant fractions

These, and other factors, make calling somatic variants difficult and still an area that is being heavily researched. If someone tells you that somatic variant calling is a solved problem, they probably have never tried to call somatic variants.

Sounds like somatic / tumor variant calling is something that will be solved by improvements at the wet lab side ( single cell selection / amplification / sequencing ) . Rather than at the computational side.

Well, single cell has a role to play (and would have more of one if WGA wasn't so lossy), but realistically, you can't sequence billions of cells from a tumor individually. Bulk sequencing still is going to have a role for quite a while.

Hell germ line calling isn't even a solved problem. Still get lots of false positives (and false negatives). It just tends to work so well that it is hard to improve it much except by making it faster, less memory intensive, etc

Solved was the wrong word. I just meant improved. There is only so much you can do at the computational side. Wet lab also has its part to play.

A germline variant caller generally has a ploidy-based genotyping algorithm built in to part of the algorithm/pipeline. I believe, IIRC, the GATK UnifiedGenotyper for instance does both variant calling and then genotype calling. So to call a genotype for a variant it is expecting a certain number of reads to support the alternative allele. When working with somatic variants all of the assumptions about how many reads you expect with a variant at a position to distinguish between true and false positives are no longer valid. Except for fixed mutations throughout the tumor population only some proportion of cells will hold a somatic variation. You also typically have some contamination from normal non-cancerous cells. Add in complications from significant genomic instability with lots of copy number variations and such and you have a need for a major change in your model for calling variation while minimizing artifactual calls. So you have a host of other programs that have been developed specifically for looking at somatic variation in tumor samples.

一篇文献:

Comparison of somatic mutation calling methods in amplicon and whole exome sequence data

是qiagen公司发的

High-throughput sequencing is rapidly becoming common practice in clinical diagnosis and cancer research. Many algorithms have been developed for somatic single nucleotide variant (SNV) detection in matched tumor-normal DNA sequencing. Although numerous studies have compared the performance of various algorithms on exome data, there has not yet been a systematic evaluation using PCR-enriched amplicon data with a range of variant allele fractions. The recently developed gold standard variant set for the reference individual NA12878 by the NIST-led “Genome in a Bottle” Consortium (NIST-GIAB) provides a good resource to evaluate admixtures with various SNV fractions.

Using the NIST-GIAB gold standard, we compared the performance of five popular somatic SNV calling algorithms (GATK UnifiedGenotyper followed by simple subtraction, MuTect, Strelka, SomaticSniper and VarScan2) for matched tumor-normal amplicon and exome sequencing data.

Nevertheless, detecting somatic mutations is still challenging, especially for low-allelic-fraction variants caused by tumor heterogeneity, copy number alteration, and sample degradation

We used QIAGEN’s GeneRead DNAseq Comprehensive Cancer Gene Panel (CCP, Version 1) for enrichment and library construction in triplicate。

QIAGEN’s GeneRead DNAseq Comprehensive Cancer Gene Panel (Version 1) was used to amplify the target region of interest (124 genes, 800 Kb).

When analyzing different types of data, use of different algorithms may be appropriate.

DNA samples of NA12878 and NA19129 were purchased from Coriell Institute. Sample mixtures were created based on the actual amplifiable DNA in each sample, resulting in 0%, 8%, 16%, 36%, and 100% of NA12878 sample mixed in the NA19129 sample, respectively.We treated the mixed samples at 8%, 16%, 36%, and 100% as the virtual tumor samples and the 0% as the virtual normal sample.

五个软件的算法是:

1. NaiveSubtract — SNVs were called separately from virtual tumor and normal samples using GATK UnifiedGenotyper [22]. For exome sequencing data, reads were already mapped, locally realigned and recalibrated by the 1,000 Genomes Project. So SNVs were directly called on the BAM files using GATK Unified Genotyper. Then, SNVs detected in the virtual normal sample were removed from the list of SNVs detected in the virtual tumor sample, leaving the “somatic” SNVs.

2. MuTect — MuTect is a method developed for detecting the most likely somatic point mutations in NGS data using a Bayesian classifier approach. The method includes pre-processing aligned reads separately in tumor and normal samples and post-processing resulting variants by applying an additional set of filters. We ran MuTect under the High-Confidence mode with its default parameter settings. We disabled the “Clustered position” filter and the “dbSNP filter” for the amplicon sequencing reads, and we disabled the “dbSNP filter” for the exome sequencing.

3. SomaticSniper — SomaticSniper calculates the Bayesian posterior probability of each possible joint genotype across the normal and cancer samples. We tuned the software’s parameters to increase sensitivity and then filtered raw results using a Somatic Score cut-off of 20 to improve specificity.

4. Strelka — Strelka reports the most likely genotype for tumor and normal samples based on a Bayesian probability model. Post-calling filters built into the software are based on factors such as read depth, mismatches, and overlap with indels. We skipped depth filtration for exome and amplicon sequencing data as recommended by the Strelka authors. For the amplicon sequencing reads, we set the minimum MAPQ score at 17 for consistency with the defaults in GATK UnifiedGenotyper. We used variants passing Strelka post-calling filters for analysis.

5. VarScan2 — VarScan2 performs analyses independently on pileup files from the tumor and normal samples to heuristically call a genotype at positions achieving certain thresholds of coverage and quality. Then, sites of the genotypes not matched in tumor and normal samples are classified into somatic, germline, or ambiguous groups using Fisher’s exact test. We generated the pileup files using SAMtools mpileup command.

The compatibility of the output VCF files between different methods as well as the NIST-GIAB gold standard was examined using bcbio.variation tools and manual inspection. The reported SNP call representations between files are comparable to each other.

来自于文献:http://www.biomedcentral.com/1471-2164/15/244

05

使用oncotator做突变注释

功能:vcf格式突变数据进一步注释成maf格式

做过癌症数据分析的童鞋都知道,TCGA里面用maf格式来记录突变!那么maf格式的数据是如何得来的呢,我们都知道,做完snp-calling一般是得到vcf格式的突变记录数据文件,然后再用annovar或者其它蛋白结构功能影响预测软件注释一下,还远达不到maf的近100条记录。

而大名鼎鼎的broad institute就规定了maf格式的突变注释文件,他就是利用了十几个常见的已知数据库来注释我们得到的vcf突变记录,通常是对somatic的突变才注释成maf格式的数据!
大名鼎鼎的broadinstitute出品的突变注释工具:http://www.ncbi.nlm.nih.gov/pubmed/25703262
本身也是一个在线工具:
集成了下面所有的分析资源
而且还提供了API

Genomic Annotations

  • Gene, transcript, and functional consequence annotations using GENCODE for hg19.
  • Reference sequence around a variant.
  • GC content around a variant.
  • Human DNA Repair Gene annotations from Wood et al.

Protein Annotations

  • Site-specific protein annotations from UniProt.
  • Functional impact predictions from dbNSFP.

Cancer Variant Annotations

Non-Cancer Variant Annotations

因为要下载的数据有点多,我这里就不用自己的电脑测试了,安装过程也很简单的!