06

拷贝数变异检测芯片介绍

这里的拷贝数变异检测芯片指的是Affymetrix Genome-Wide Human SNP Array 6.0

cel数据,需要处理成segment及genotype数据
这个芯片在TCGA计划里面用的非常多,是标配了。大家只要记住,这是一个跟拷贝数变异检测相关的芯片,而且还可以测一些genotype  
Affymetrix Genome-Wide Human SNP Array 6.0是唯一可以真正将CNP(拷贝数多态性)转化成高分辨率的参考图谱的平台。主要应用领域包括全基因组SNP分型、全基因组CNV分型、全基因组关联 分析、全基因组连锁分析。除了进行基因分型外,还为拷贝数研究和LOH研究提供帮助,从而能够进行:UPD检测、亲子鉴定、异常的亲代起源分析(针对 UPD和缺失)、纯合性分析、血缘关系鉴定。
SNP Array 6.0是昂飞公司继Mapping10k、100k、500k和SNP5.0芯片后推出的新一代SNP芯片。在一张芯片上可以分析一个样本906,600 个SNP的基因型, 大约有482,000个SNP来自于前代产品500K和SNP5.0芯片。剩下424,000个SNP包括了来源于国际HapMap计划中的标签 SNP,X,Y染色体和线粒体上更具代表性的SNP,以及来自于重组热点区域和500K芯片设计完成后新加入dbSNP数据库的SNP。该芯片同时含 946,000个非多态性CNV探针,用于检测拷贝数变异,其中202,000个用于检测5677个已知拷贝数变异区域的探针,这些区域来源于多伦多基因 组变异体数据库。该数据库中每隔3,182个非重叠片段区域分别用61个探针来检测。除了检测这些已知的拷贝数多态区域,还有超过744,000个探针平 均分配到整个基因组上,用来发现未知的拷贝数变异区域。SNP和CNV两种探针高密度且均匀地分布在整个基因组,作为拷贝数变异和杂合性缺失(LOH)检 测的工具来发现微小的染色体增加和缺失。为广大生命科学研究者提高发现复杂疾病相关基因的可能提供了强有力的工具。
通过与哈佛大学合办的Broad研究所合作,SNP6.0芯片在数据准确性和一致性方面达到了新的高度。相应推出的Genotyping Console用来处理SNP6.0芯片数据和全基因组遗传分析及质量控制。

产品特点:

1.涵盖超过1,800,000个遗传变异标志物:包括超过906,600个SNP和超过946,000个用于检测拷贝数变化(CNV,Copy Number Variation)的探针;

2.SNP和CNV两种探针高密度且均匀地分布在整个基因组,不仅可以用于SNP基因精确分型,还可用于拷贝数变异CNV的研究;

3.744,000个探针平均分配到整个基因组上,用来发现未知的拷贝数变异区域;

4.可用于Copy-neutral LOH/UPD检测,亲子鉴定,纯合性分析、血缘关系鉴定、遗传病或其它疾病的研究。

参考:http://www.biomart.cn/specials/cnv2014/article/84169

在NCBI的GEO数据库里面可以查到这个芯片,已经有一万多个样本数据啦!
图中第一个是CCLE计划的近千个样本,可能是定制化了的snp6.0芯片吧
clipboard
使用这个芯片数据来发文章的非常多,见列表:http://media.affymetrix.com/support/technical/other/snp6_array_publications.pdf
还有一篇2010-nature文章讲了如何用picnic来研究cnv,http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3145113/
也有一篇2010年的文章提出了新的软件来分析这个芯片cnv数据http://bioinformatics.oxfordjournals.org/content/26/11/1395.long
实现同样功能的软件,非常之多,还有一个R的bioconductor系列的包
clipboard2
随便进去都可以找到很多raw data,可以自己进行分析的!

 

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

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

 

十二 30

用GSEA来做基因集富集分析

how to use GSEA?
这个有点类似于pathway(GO,KEGG等)的富集分析,区别在于gene set(矫正好的基于文献的数据库)的概念更广泛一点,包括了

how to download GSEA ?

what's the input for the GSEA?

说明书上写的输入数据是:GSEA supported data files are simply tab delimited ASCII text files, which have special file extensions that identify them. For example, expression data usually has the extension *.gct, phenotypes *.cls, gene sets *.gmt, and chip annotations *.chip. Click the More on file formats help button to view detailed descriptions of all the data file formats.
实际上没那么复杂,一个表达矩阵即可!然后做一个分组说明的cls文件即可。
主要是自己看说明书,做出要求的数据格式:http://www.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats
表达矩阵我这里下载GSE1009数据集做测试吧!
cls的样本说明文件,就随便搞一搞吧,下面这个是例子:
6 2 1
# good bad
good good good bad bad bad
文件如下,六个样本,根据探针来的表达数据,分组前后各三个一组。
clipboard
现在开始运行GSEA!

start to run the GSEA !

首先载入数据
clipboard
确定无误,就开始运行,运行需要设置一定的参数!
clipboard

what's the output ?

输出的数据非常多,对你选择的gene set数据集里面的每个set都会分析看看是否符合富集的标准,富集就出来一个报告。
点击success就能进入报告主页,里面的链接可以进入任意一个分报告。
最大的特色是提供了大量的数据集:You can browse the MSigDB from the Molecular Signatures Database page of the GSEA web site or the Browse MSigDB page of the GSEA application. To browse the MSigDB from the GSEA application:
 
有些文献是基于GSEA的:

 

十二 29

用firehose_get 来下载所有TCGA寄存在broad的数据

该软件是broad institute写的一个数据接口,主要是供他人下载TCGA的所有寄存在broad institute的免费数据,主要是level3,level4的数据。(说错了,好像只有level4的数据,就是可以发文章的分析结果及图片)
软件下载地址:https://confluence.broadinstitute.org/display/GDAC/Download

懂它的使用规则,编码规则即可:
就是一个很简单的shell脚本而已,根据几个用户自定义参数来选择性的下载数据。
clipboard
我们可以用-t这个参数来指定下载的数据类型,可以是mut/rna/mutsig/gistic等各种数据,至于这些单词代表什么意义,需要自己去看说明书啦
还可以指定时间,截止到什么时间的数据!
它支持的癌症种类:

ACC  BLCA  BRCA  CESC  COAD  COADREAD  DLBC  ESCA  
	GBM  HNSC  KICH  KIRC  KIRP  LAML  LGG  LIHC  
	LUAD  LUSC  OV  PAAD  PANCANCER  PANCAN8  PANCAN12  PRAD  
	READ  SARC  SKCM  STAD  THCA  UCEC  UCS
这些癌症种类的简称,也是可以去官网里面看到的!官网:http://gdac.broadinstitute.org

 

十二 29

所有TCGA收集的mRNA表达数据集数据集-GSE62944

无意中看一篇文献,有提到这个数据集,我简单看了一下,居然还真的这么全面!!!
它处理了目前(大概是2015年6月)TCGA收集的所有癌症样本的mRNA表达数据,并且统一处理成了count和RPKM两种表达量形式。
GEO地址:http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62944

Title Alternatively processed and compiled RNA-Sequencing and clinical data for thousands of samples from The Cancer Genome Atlas
Organism Homo sapiens
Experiment type Expression profiling by high throughput sequencing
Summary We reprocessed RNA-Seq data for 9264 tumor samples and 741 normal samples across 24 cancer types from The Cancer Genome Atlas with "Rsubread". Rsubread is an open source R package that has shown high concordance with other existing methods of alignment and summarization, but is simple to use and takes significantly less time to process data. Additionally, we provide clinical variables publicly available as of May 20, 2015 for the tumor samples where the TCGA ids are matched.
这样就非常方便大家使用了。
可以直接拿来进行聚类,看看不同癌症种类如何聚集在一起,也可以直接拿来跟某个临床指标进行关联分析.
比如,如果我们想用DEseq来做差异分析,那么需要的数据就是基因的count,这里就有近一万个癌症样本的基因表达的count值,随便都可以分类看看差异情况,再富集分析分析。非常适合初学者练手。
癌症种类见官网:http://gdac.broadinstitute.org/
clipboard

 

十二 29

自动化出网页报告的R语言包-Nozzle

根据测试代码输出的报告如下:
clipboard
这些报告里面的元素,都是测试代码添加的,很容易就能理解
十二 29

用Mutation-Assessor软件来看突变位点对基因或者蛋白功能的影响

这是一个在线工具,非常好用,对snp位点进行注释,看看该突变是否影响蛋白功能,一定要收藏!!!

官网:http://mutationassessor.org/

也应该是有standalone版本,我没有去找,不过网页版就很好用了,只需要复制粘贴进去自己想分析的数据,按照一定的格式即可,比如:
clipboard
该软件就能分析出该突变位点发生在BRCA2这个基因上面,对氨基酸的改变也能写出来,对蛋白功能的改变等选项都是可以自由定制化的。
输入数据非常广泛:
The server accepts list of variants, one variant per line, plus optional text describing your variants,
in genomic coordinates, "+" strand assumed :
<genome build>,<chromosome>,<position>,<reference allele>,<substituted allele>
Genome build is optional (build 18 assumed), accepted values: 'hg18' and 'hg19'
Examples:

hg19,13,32912555,G,T   BRCA2
hg18,7,55178574,G,A   GBM
7,55178574,G,A   GBM

or in protein space: <protein ID> <variant> <text>, where <protein ID> can be :

1. Uniprot protein accession (i.e. EGFR_HUMAN)
2. NCBI Refseq protein ID (i.e. NP_005219)

examples:

EGFR_HUMAN R521K
EGFR_HUMAN R98Q Polymorphism
EGFR_HUMAN G719D disease
NP_000537 G356A
NP_000537 G360A dbSNP:rs35993958
NP_000537 S46A Abolishes phosphorylation

ID types can be mixed in one list in any way.

 

十二 25

使用进度条和并行计算来模拟计算赌博输赢概率

http://www.zhihu.com/question/19605599

以前看到过一个赌徒「必胜方法」
有一个广泛流传于赌徒中的“必胜方法”,是关于赌徒悖论极好的说明:一个赌场设有“猜大小”的赌局,玩家下注后猜“大”或者猜“小”,如果输了,则失去赌注,如果赢的话,则获得本金以及0.9倍的利润。 “必胜法”是这么玩的:
(1)押100猜“大”;
(2)如果赢的话,返回(1)继续;
(3)如果输的话则将赌注翻倍后继续猜“大”,因为不可能连续出现“大”,总会有获胜的时候,而且由于赌注一直是翻倍,只要赢一次,就会把所有输掉的钱赢回;
(4)只要赢了,就继续返回(1)。

一个朋友说写了程序来验证正确与否,我感觉蛮有趣的,想想自己也学了不少,也应该实践一下啦!
写出程序不难,就一个递归而已,不过是应该用并行计算,不然算的太慢了
另外一个重点是,如何随机???
[R]
judge=function(l=left,c=count,n=number,w=win){
c=c+1;
if(rnorm(1)>0){
l=l+2^n
n=1 # win
w=w+1
}else{
l=l-2^n
n=n+1 #lose
}
return(c(l,c,n,w))
}
my_fun=function(x){
tmp=judge(0,1,1,0)
for(i in 1:x){
tmp=judge(tmp[1],tmp[2],tmp[3],tmp[4])
print(tmp)
}
return(tmp[1])
}
my_fun(10)
[/R]
第一个函数是根据现在的状态,接受你的本金,第几次赌了,以及第几次连续下注,最后,你总共赢了多少次,这四个参数,然后随机给一个大小概率,这样你的本金会变化,赌的次数会加1
第二个函数是传入我们需要赌多少次,看结果:
> my_fun(10)
[1] -2  2  2  0
[1] 2 3 1 1
[1] 0 4 2 1
[1] -4  5  3  1
[1] 4 6 1 2
[1] 6 7 1 3
[1] 8 8 1 4
[1] 6 9 2 4
[1]  2 10  3  4
[1] 10 11  1  5
[1]  8 12  2  5
[1] 8
如果我们赌11次,可以看到,我最后会剩余8块钱,每次输赢的情况都反应在里面了,可以自己模拟多次看看!
因为我只赌了11次,所以很快,如果我赌1000次,而且还想检验一下10000次模拟结果,就会比较慢了!
我首先使用进度条模拟一下结果,代码如下:
2
还是比较慢的##Time difference of 1.861802 mins
我用了apply,好像时间是节省了一些,不过聊胜于无!
> library(pbapply)
> start=Sys.time()
> results=pbsapply(1:10000,function(i){my_fun(1000)})
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
> Sys.time()-start  ##9.909524 secs
Time difference of 1.301539 mins
那接下来我们分析一下模拟的结果吧
> summary(results)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -64580     974     998     985    1020    1110
可以看到,我们平均下来是会赢钱的,而且赢面很大,赢的次数非常多!!!
但是,我们看下面的图就知道,一个很诡异的结果!
而且,我这里用的模拟胜负,并不是很好
3
这其实还只是小批量模拟,如果要模拟百亿次,首先我的笔记本肯定不行,cpu太破了,如果用服务器,就需要用并行计算啦!
##下面是用多核并行计算的代码,大家有兴趣可以自己去玩一下
library(parallel)
cl.cores <- detectCores() 
cl <- makeCluster(cl.cores)
clusterExport(cl, "judge") 
start=Sys.time()
results=parSapply( 
  cl=cl, 
  1:10000,
  my_fun(1000)
)
Sys.time()-start  ##4.260994 secs
stopCluster(cl)


十二 25

做癌症研究一定要把这几十篇TCGA的大文章看完

都是发表在nature,cell还有新英格兰医学杂志上面的超级文章!每个文章附件都有一百多页,比博士论文还长,但是它们的分析套路其实都一样,都是那几种数据,包括WGS,WES,RNA-Seq,芯片表达量,miRNA表达量,甲基化数据,蛋白数据。分析过程也差不多,无法就是对癌症进行进一步的分类,癌症亚型,或者看看driver mutation,进一步解释癌症病变,转移,扩散机理,或者找标记物signature,辅助治疗等等,具体的要等我把这些文献看完了才能再进一步讲解,请做癌症研究方向的一定要把它们看完。

1

我已经下载完了,大家如果没有权限下载,就需要自己想办法啦!

image

非常值得大家阅读!!!

 

十二 24

使用R包cgdsr来下载TCGA的数据

前面我讲到TCGA的数据可以在5个组织机构可以获取,他们都提供了类似的接口来供用户下载数据

每个接口都有使用教程,比如http://firebrowse.org/tutorial/FireBrowse-Tutorial.pdf

非常详细!!!

有的还专门写了软件接口:https://confluence.broadinstitute.org/display/GDAC/Download

或者写了R的接口:http://www.cbioportal.org/cgds_r.jsp

接下来我们要讲的就是cbioportal网站提供的一个R接口,非常好用,只需记住4个函数即可!!! Continue reading

十二 24

TCGA数据下载大全

并不是所有的数据都能下载,很多数据需要有权限才能下载的!!!
首先我们可以根据TCGA的文章来下载数据:

总共也就几十篇文章,都是发表在大杂志上面的。
每篇文章都会提供数据的打包下载,例如:

The molecular taxonomy of primary prostate cancer
Cell Volume 163 Issue 4: p1011-1025 Read the full article
Portal Publication Site and Associated Data Files

Comprehensive Molecular Characterization of Papillary Renal Cell Carcinoma
NEJM. Published on line on Nov 4th, 2015 Read the full article
Portal Publication Site and Associated Data Files

那个portal链接点击进去,就可以看到所有的下载链接了!
这是根据文章来分类打包好的数据!

同时也可以通过其它数据接口来下载

Tools for Exploring Data and Analyses

TCGA Data Portal

这几个接口都挺好用的:
非常详细!!!而且还专门写了软件接口:https://confluence.broadinstitute.org/display/GDAC/Download
或者写了R的接口:http://www.cbioportal.org/cgds_r.jsp
一般都推荐用TCGA自己的数据接口:https://tcga-data.nci.nih.gov/tcga/
里面对所有的样本都进行了统计
通过https://tcga-data.nci.nih.gov/tcga/dataAccessMatrix.htm可以进行定制化的数据下载!
里面有很多TCGA自定义的名词:

Data Levels and Data Types

https://tcga-data.nci.nih.gov/docs/dictionary/ 可以看到所有名词的解释:
数据的种类如下:
还记得以前看到一篇TCGA自己的关于胃癌的文章,发表在nature上面,文章涉及到了TCGA的各个方面的分析,所以附件PDF竟然有133页!!!

包含的其它数据有:
十二 23

看看Y染色体上面的基因在测序数据里的覆盖度和测序深度

作用:可以检测别人是否把自己的样本搞混,也可以看看测序是否分布均匀!
首先,我们要拿到Y染色体上面的基因的坐标信息!
因为我们的是hg19,所以我们要下载hg19的基因信息!
我们首先解析refGene文件,找到chrY的unique基因!
这四列分别是:chromosome/start/end/gene_symbol
clipboard
4程序如下:

[perl]
open FH,"/home/jmzeng/hg19/chrY.gene.special.position" or die "file error !!!";
while(<FH>){
    chomp;
    @F=split;
    foreach ($F[1]..$F[2]){
        $h{$_}=$F[3];
    }
    $length{$F[3]}=$F[2]-$F[1]+1;
}
close FH;
open FH,$ARGV[0];
while(<FH>){
    chomp;
    @F=split;
    next unless $F[0] eq 'chrY';
    next if $F[2]<20;
    if (exists $h{$F[1]}){
        $count{$h{$F[1]}}++ ;
    }else{
        $count{'other'}++ ;
    }    
}
close FH;
print "$_\t$length{$_}\t$count{$_}\n" foreach sort keys %count;</pre>
</div>
<div>[/perl]

对一个男性样本,结果会如下:
gene/length/pos
AMELY 8111 1269
BCORP1 47724 689
CSPG4P1Y 3799 538
DAZ1 69739 762
DAZ2 71901 228
DAZ3 73222 233
DAZ4 73222 540
DDX3Y 12825 3654
EIF1AY 17445 929
FAM224A 4295 82
FAM224B 4293 85
GOLGA2P3Y 4866 68
GYG2P1 15476 547
HSFY2 42277 3950
KDM5D 39526 7425
NLGN4Y 319396 3872
PCDH11Y 105374 6627
PRKY 107577 1390
PRORY 3388 735
RBMY1B 14451 232
RBMY1D 14411 117
RBMY1E 14410 157
RBMY1J 14407 65
RBMY2EP 6416 27
RBMY2FP 7348 419
RPS4Y1 25376 1856
RPS4Y2 24966 1831
SRY 888 703
TBL1Y 180999 3231
TGIF2LY 958 808
TMSB4Y 2457 534
TSPY4 132211 1525
TTTY14 205048 394
TTTY4C 36811 39
TTTY9A 9317 580
TXLNGY 23067 1968
USP9Y 159610 10508
UTY 232293 6670
VCY 742 291
XKRY2 1582 980
ZFY 47437 3125
other 100328
对女性样本,结果会如下;
NLGN4Y 319396 575
PCDH11Y 105374 1643
PRKY 107577 82
TGIF2LY 958 191
TTTY14 205048 139
other 54297
从结果可以看出来,很多基因都是y染色体特有的,这个结果是表明我们的测序非常棒

 

十二 22

根据chrY独有区域的覆盖度及测序深度来判断性别

这个也是基于bam文件来的,判断chrY独有区域的覆盖度及测序深度
首先下载chrY独有区域的记录文件,https://www.familytreedna.com/documents/bigy_targets.txt
然后用samtools depth来统计测序深度,samtools  depth $i |grep 'chr[XY]'
depth统计结果文件如下:
mzeng@ubuntu:/home/jmzeng/gender_determination$ head Sample3.depth
chrX    60085    1
chrX    60086    1
chrX    60087    1
chrX    60088    1
chrX    60089    1
chrX    60090    1
chrX    60091    1
chrX    60092    1
chrX    60093    1
chrX    60094    1
然后我随便写了一个脚本来对测序深度文件进行再统计,统计覆盖度及测序深度

[perl]
open FH,"bigy_targets.txt";
while(<FH>){
 chomp;
 @F=split;
 $all+=$F[2]-$F[1]+1;
 foreach ($F[1]..$F[2]){
  $h{$_}=1;
 }
}
close FH;
open FH,$ARGV[0];
while(<FH>){
 chomp;
 @F=split;
 next unless $F[0] eq 'chrY';
 if (exists $h{$F[1]}){
  $pos++;
  $depth+=$F[2];
 }
}
close FH;
$average=$depth/$pos;
$coverage=$pos/$all;
print "$pos\t$average\t$coverage\n" ;

[/perl]

 

 
这样对那三个样本结果如下:
clipboard
可以看到只有sample4,是覆盖率极低的,而且记录到的pos位点也特别少,所以她是女性!
这里测序深度没有意义。
十二 22

根据X,Y染色体比对上的reads数来判断性别

针对高通量测序数据,包括WGS,WES
我这里主要讲根据bam文件里面的chrX和chrY的reas比例来判断性别,大家可以自己处理数据得到bam文件!
主要是读取bam文件,选择chrX上面的记录,统计genotype即可:
也可以统计测序深度,samtools  depth $i |grep 'chr[XY]'
如果chrX和chrY的reads比例超过一定值,比如50倍,就判定为女性!
脚本很简单,就统计bam文件的第三列就可以啦,第三列就是染色体信息
 samtools  view $i |  perl -alne '{$h{$F[2]}++}END{print "$_\t$h{$_}" foreach sort keys %h}' >$out.chromosome.stat
对于Sample3,
chrX 1233061
chrY 140506
 
对于Sample4,很明显,这个是女性!
chrX 2734815
chrY 51860
 
对于Sample5
chrX 1329302
chrY 156663
十二 22

根据X染色体的snp的纯和率来判断性别

针对高通量测序数据,包括WGS,WES,甚至snp6芯片也行。
我这里主要讲根据vcf文件里面的chrX的纯合率来判断性别,大家可以自己处理数据得到vcf文件!
主要是读取vcf文件,选择chrX上面的记录,统计genotype即可:
clipboard
我这里拿之前的自闭症项目数据来举例子:
根据数据提供者的信息,3-4-5分别就是孩子、父亲、母亲,统计chrX的snp的的纯合和杂合的比例,代码很简单
vcf文件一般第一列是染色体,第6列是质量,第10列是基因型已经测序深度相关信息

cat Sample5.gatk.UG.vcf |perl -alne '{next unless $F[0] eq "chrX";next unless $F[5]>30;$h{(split/:/,$F[9])[0]}++}END{print "$_\t$h{$_}" foreach keys %h}' 

如果纯合的snp是杂合的倍数超过一定阈值,比如4倍,就可以判断是男性。

对于Sample3来说,很明显,是男孩,因为X染色体都是纯合突变

0/1 391
1/1 2463
2/2 1
1/2 32
0/2 1
对于Sample4来说,很明显,应该是母亲,证明之前别人给我的信息有误
1/1 3559
1/2 27
0/1 1835
0/2 5
那么Sample5很明显就是父亲咯
1/1 2626
0/1 356
1/2 22
十二 21

用重抽样+主成分方法来做富集分析

之前我们用超几何分布检验的方法做了富集分析,使用的是GSE63067.diffexp.NASH-normal.txt的logFC的绝对值大于0.5,并且P-value小雨0.05的基因作为差异基因来检验kegg的pathway的富集情况

结果是这样的

image001

我们接下来用另外一种方法来做富集分析,顺便检验一下,是不是超几何分布统计检验的富集分析方法就是最好的呢?

这种方法是-重抽样+主成分分析

大概的原理是,比如对上图中的,04380这条pathway来说,总共有128个基因,那么我从原来的表达矩阵里面随机抽取128个基因的表达矩阵做主成分分析,并且抽取一千次,每次主成分分析都可以得到第一主成分的贡献度值。那么,当我并不是随机抽取的时候,我就抽04380这条pathway的128个基因,也做主成分分析,并且计算得到第一主成分分析的重要性值。我们看看这个值,跟随机抽1000次得到的值差别大不大。

这时候就需要用到表达矩阵啦!

setwd("D:\\my_tutorial\\补\\用limma包对芯片数据做差异分析")

exprSet=read.table("GSE63067_series_matrix.txt.gz",comment.char = "!",stringsAsFactors=F,header=T)

rownames(exprSet)=exprSet[,1]

exprSet=exprSet[,-1]

我们根据ncbi里面对GSE63067的介绍可以知道,对应NASH和normal的样本的ID号,就可以提取我们需要的表达矩阵

把前面两属于Steatosis的样本去掉即可,exprSet=exprSet[,-c(1:2)]

然后再把芯片探针的id转换成entrez id

exprSet=exprSet[,-c(1:2)]

library(hgu133plus2.db)

library(annotate)

platformDB="hgu133plus2.db";

probeset <- rownames(exprSet)

rowMeans <- rowMeans(exprSet)

EGID <- as.numeric(lookUp(probeset, platformDB, "ENTREZID"))

match_row=aggregate(rowMeans,by=list(EGID),max)

colnames(match_row)=c("EGID","rowMeans")

dat=data.frame(EGID,rowMeans,probeset)

tmp_prob=merge(dat,match_row,by=c("EGID","rowMeans"))

relevantProbesets=as.character(tmp_prob$probeset)

length(relevantProbesets) #hgu133plus2.db  20156

exprSet=exprSet[relevantProbesets,]

EGID_name=as.numeric(lookUp(relevantProbesets, platformDB, "ENTREZID"))

rownames(exprSet)=as.character(EGID_name)

d=exprSet

最后得到表达矩阵表格

image002

我们首先得到1000次随机挑选128个基因的表达矩阵的主成分分析,第一主成分贡献度值。

gene128=sapply(1:1000,function(y) {

dat=t(d[sample(row.names(d), 128, replace=TRUE), ]);

round(100*summary(fast.prcomp(dat))$importance[2,1],2)

}

)

很快就能得到结果,可以看到数据如下

>  summary(gene128)

Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

19.1    25.8    28.8    29.8    32.5    59.7

image003

 

那么接下来我们挑选这个04380这条pathway特有128个基因来算第一主成分贡献度值

path_04380_gene=intersect(rownames(d),as.character(Path2GeneID[['04380']]))

dat=t(d[path_04380_gene,]);

round(100*summary(fast.prcomp(dat))$importance[2,1],2)

image004

得到的值是38.83,然后看看我们的这个38.83在之前随机得到的1000个数里面是否正常,就按照正态分布检验来算

1-pnorm((38.83-mean(gene128))/sd(gene128))

[1] 0.0625

可以看到已经非常显著的不正常了,可以说明这条通路被富集了。

至少说明超几何分布检验的方法得到的富集分析结果跟我们这次的重抽样+主成分分析结果是一致的,当然,也有不一致的,不然就不用发明一种新的方法了。

如果写一个循环同样可以检验所有的通路,但是这样就不需要事先准备好差异基因啦!!!这是这个分析方法的特点!

十二 21

主成分分析略讲

主成分分析是为了简化变量的个数。
我这里不涉及到任何高级统计知识来简单讲解一下主成分分析,首先我们用下面的代码随机创造一个矩阵:

options(digits = 2)
x=c(rnorm(5),rnorm(5)+4)
y=3*c(rnorm(5),rnorm(5)+4)
dat=rbind(x,y,a=0.1*x,b=0.2*x,c=0.3*x,o=0.1*y,p=0.2*y,q=0.3*y)
colnames(dat)=paste('s',1:10,sep="")
dat
library(gmodels)
pca=fast.prcomp(t(dat))
pca
summary(pca)$importance
biplot(pca, cex=c(1.3, 1.2));

Continue reading

十二 15

用超几何分布检验做富集分析

我们可以直接使用R的bioconductor里面的一个包,GOstats里面的函数来做超几何分布检验,看看每条pathway是否会富集

我们直接读取用limma包做好的差异分析结果

setwd("D:\\my_tutorial\\补\\用limma包对芯片数据做差异分析")

DEG=read.table("GSE63067.diffexp.NASH-normal.txt",stringsAsFactors = F)

View(DEG)

image001

我们挑选logFC的绝对值大于0.5,并且P-value小雨0.05的基因作为差异基因,并且转换成entrezID

probeset=rownames(DEG[abs(DEG[,1])>0.5 & DEG[,4]<0.05,])

library(hgu133plus2.db)

library(annotate)

platformDB="hgu133plus2.db";

EGID <- as.numeric(lookUp(probeset, platformDB, "ENTREZID"))

length(unique(EGID))

#[1] 775

diff_gene_list <- unique(EGID)

这样我们的到来775个差异基因的一个list

首先我们直接使用R的bioconductor里面的一个包,GOstats里面的函数来做超几何分布检验,看看每条pathway是否会富集

library(GOstats)

library(org.Hs.eg.db)

#then do kegg pathway enrichment !

hyperG.params = new("KEGGHyperGParams", geneIds=diff_gene_list, universeGeneIds=NULL, annotation="org.Hs.eg.db",

categoryName="KEGG", pvalueCutoff=1, testDirection = "over")

KEGG.hyperG.results = hyperGTest(hyperG.params);

htmlReport(KEGG.hyperG.results, file="kegg.enrichment.html", summary.args=list("htmlLinks"=TRUE))

结果如下:

image002

但是这样我们就忽略了其中的原理,我们不知道这些数据是如何算出来的,只是由别人写好的包得到了结果罢了。

事实上,这个包的这个hyperGTest函数无法就是包装了一个超几何分布检验而已。

如果我们了解了其中的统计学原理,我们完全可以写成一个自建的函数来实现同样的功能。

超几何分布很简单,球分成黑白两色,数量已知,那么你随机抽有限个球,应该抽多少白球的问题!
公式就是 exp_count=n*M/N
然后你实际上抽了多少白球,就可以计算一个概率值!
换算成通路的富集概念就是,总共有多少基因,你的通路有多少基因,你的通路被抽中了多少基因(在差异基因里面属于你的通路的基因),这样的数据就足够算出上面表格里面所有的数据啦!
tmp=toTable(org.Hs.egPATH)
GeneID2Path=tapply(tmp[,2],as.factor(tmp[,1]),function(x) x)
Path2GeneID=tapply(tmp[,1],as.factor(tmp[,2]),function(x) x)
#phyper(k-1,M, N-M, n, lower.tail=F)
#n*M/N
diff_gene_has_path=intersect(diff_gene_list,names(GeneID2Path))
n=length(diff_gene_has_path) #321 # 这里算出你总共抽取了多少个球
N=length(GeneID2Path) #5870  ##这里算出你总共有多少个球(这里是错的,有多少个球取决于背景基因!一般是两万个)
options(digits = 4)
for (i in names(Path2GeneID)){
 M=length(Path2GeneID[[i]])  ##这个算出你的所有的球里面,白球有多少个
 exp_count=n*M/N  ###这里算出你抽取的球里面应该多多少个是白色
 k=0         ##这个k是你实际上抽取了多少个白球
 for (j in diff_gene_has_path){
 if (i %in% GeneID2Path[[j]]) k=k+1
 }
 OddsRatio=k/exp_count
 p=phyper(k-1,M, N-M, n, lower.tail=F)  ##根据你实际上抽取的白球个数,就能算出富集概率啦!
 print(paste(i,p,OddsRatio,exp_count,k,M,sep="    "))
}
随便检查一下,就知道结果是一模一样的!

 

十二 15

下载所有的酶的信息,并且解析好

之前我提到过kegg数据库里面的有些pathway下面没有对应任何基因,当时我还在奇怪,怎么会有这样的通路呢?

然后,我随机挑选了其中一条通路(hsa01000),进行查看,发现正好是所有的酶的信息。

好奇怪,不明白为什么kegg要列出所有的酶信息

http://www.genome.jp/kegg-bin/get_htext?hsa01000

image001

下载htext格式的酶的信息

798K Dec 15  2015 hsa01000.keg

查看文件,发现也是层级非常清楚的结构,D已经是最底级别的酶了,而E对应的基因是属于该酶的。

简单统计一下,发现跟酶相关的基因有3688个,而最底级别的酶有5811个,应该会持续更新的

image002

如果想做成kegg那样的基因与酶对应表格,也是非常简单的!

 

 

国际系统分类法按酶促反应类型,将酶分成六个大类:

1、氧化还原酶类(oxidoreductases) 催化底物进行氧化还原反应的酶类。包括电子或氢的转移以及分子氧参加的反应。常见的有脱氢酶、氧化酶、还原酶和过氧化物酶等。

2、转移酶类(transferases) 催化底物进行某些基团转移或交换的酶类,如甲基转移酶、氨基转移酶、转硫酶等。

3、水解酶类(hydrolases)催化底物进行水解反应的酶类。如淀粉酶、粮糖苷酶、蛋白酶等。

4、裂解酶类(lyases)或裂合酶类(synthases) 催化底物通过非水解途径移去一个基团形成双键或其逆反应的酶类,如脱水酶、脱羧酸酶、醛缩酶等。如果催化底物进行逆反应,使其中一底物失去双键,两底物间形成新的化学键,此时为裂合酶类。

5、异构酶类(isomerases)催化各种同分异构体、几何异构体或光学异构体间相互转换的酶类。如异构酶、消旋酶等。

6、连接酶类(ligases)或合成酶类(synthetases)催化两分子底物连接成一个分子化合物的酶类。

上述六大类酶用EC(enzyme commission)加1.2.3.4.5.6编号表示,再按酶所催化的化学键和参加反应的基团,将酶大类再进一步分成亚类和亚-亚类,最后为该酶在这亚-亚类中的排序。如α淀粉酶的国际系统分类纩号为:EC3.2.1.1

 

EC3——Hydrolases 水解酶类

EC3.2——Glycosylases 转葡糖基酶亚类

EC3.2.l——Glycosidases 糖苷酶亚亚类i.e.enzymes hydmlyzing O-and S-glycosyl compound即能水解O-和S-糖基化合物

EC3.2.1.1 Alpha-amylase, α-淀粉酶

值得注意的是,即使是同一名称和EC编号,但来自不同的物种或不同的组织和细胞的同一种酸,如来自动物胰脏、麦芽等和枯草杆菌BF7658的α-淀粉酶等,它们的一级结构或反应机制可解不同,它们虽然都能催化淀扮的水解反应,但有不同的活力和最适合反应条件。

可以按照酶在国际分类编号或其推荐名,从酶手册(Enzyme Handbook)、酶数据库中检索到该酶的结构、特性、活力测定和Km值等有用信息。著名的手册和数据库有:

手册:

1、Schomburg,M.Salzmann and D.Stephan:Enzyme Handbook 10 Volumes

2、美国Worthington Biochemical Corporation:Enzyme Manual

(http://www.worthington-biochem.com/index/manual.htm/)

数据库:

l、德国BRENDA:Enzyme Database(http://www.brenda wnzymes.org)

2、Swissprot:EXPASYENZYME Enzyme nomenclature database (http://www.expasy.org/enzyme/)

3、IntEnz:Integrated relational Enzyme database (http://www.ebi.ac.uk/mtenz)