05

自学CHIP-seq分析第六讲~寻找peaks

CHIP-seq测序的本质还是目标片段捕获测序,跟WES不同的是,它不是通过固定的芯片探针来固定的捕获基因组上面特定序列,而是根据你选择的IP不同,你细胞或者机体状态不同,捕获到的序列差异很大!而我们研究的重点,就是捕获到的差异。而我们对CHIP-seq测序数据寻找peaks的本质就是得到所有测序数据比对在全基因组之后在正个基因组上面的测序深度里面寻找比较突出的。比如对WES数据来说,各个外显子,或者外显子的5端到3端,理论上测序深度应该是一致的,都是50X~200X,画一个测序深度曲线,应该是近似于一条直线。对我们的CHIP-seq测序数据来说,在所捕获的区域上面,理论上测序深度是绝对不一样的,应该是近似于一个山峰。而那些覆盖度高的地方,山顶,就是我们的IP所结合的热点,也就是我们想要找的peaks,在IGV里面看到大致是下面这样:

chipseq-0-peakdetection-large-IGV

可以看到测序的reads分布是绝对的不均匀的!我们通常说的CHIP-seq测序的IP,可以是各个组蛋白的各个修饰位点对应的抗体,或者是各种转录因子的抗体,等等

如何定义热点呢?通俗地讲,热点是这样一些位置,这些位置多次被测得的read所覆盖(我们测的是一个细胞群体,read出现次数多,说明该位置被TF结合的几率大)。那么,read数达到多少才叫多?这就要用到统计检验喽。假设TF在基因组上的分布是没有任何规律的,那么,测序得到的read在基因组上的分布也必然是随机的,某个碱基上覆盖的read的数目应该服从二项分布。

具体统计学原理直接看原创吧:http://www.plob.org/2014/05/08/7227.html

Continue reading

05

自学CHIP-seq分析第五讲~测序数据比对

比对本质是是很简单的了,各种mapping工具层出不穷,我们一般常用的就是BWA和bowtie了,我这里就挑选bowtie2吧,反正别人已经做好了各种工具效果差异的比较,我们直接用就好了,代码如下:

## step5 : alignment to hg19/ using bowtie2 to do alignment
## ~/biosoft/bowtie/bowtie2-2.2.9/bowtie2-build ~/biosoft/bowtie/hg19_index /hg19.fa ~/biosoft/bowtie/hg19_index/hg19
## cat >run_bowtie2.sh
ls *.fastq | while read id ;
do
echo $id
#~/biosoft/bowtie/bowtie2-2.2.9/bowtie2 -p 8 -x ~/biosoft/bowtie/hg19_index/hg19 -U $id -S ${id%%.*}.sam 2>${id%%.*}.align.log;
#samtools view -bhS -q 30 ${id%%.*}.sam > ${id%%.*}.bam ## -F 1548 https://broadinstitute.github.io/picard/explain-flags.html
# -F 0x4 remove the reads that didn't match
samtools sort ${id%%.*}.bam ${id%%.*}.sort ## prefix for the output
# samtools view -bhS a.sam | samtools sort -o - ./ > a.bam
samtools index ${id%%.*}.sorted.bam
done

这个索引~/biosoft/bowtie/hg19_index/hg19需要自己提取建立好,见前文

初步比对的sam文件到底该如何过滤,我查了很多文章都没有给出个子丑寅卯,各执一词,我也没办法给大家一个标准,反正我测试了好几种,看起来call peaks的差异不大,就是得不到文章给出的那些结果!!

一般来说,初步比对的sam文件只能选取unique mapping的结果,所以我用了#samtools view -bhS -q 30,但是结果并没什么改变,有人说是peak caller这些工具本身就会做这件事,所以取决于你下游分析所选择的工具。

给大家看比对的日志吧:

SRR1042593.fastq
16902907 reads; of these:
16902907 (100.00%) were unpaired; of these:
667998 (3.95%) aligned 0 times
12467095 (73.76%) aligned exactly 1 time
3767814 (22.29%) aligned >1 times
96.05% overall alignment rate
[samopen] SAM header is present: 93 sequences.
SRR1042594.fastq
60609833 reads; of these:
60609833 (100.00%) were unpaired; of these:
9165487 (15.12%) aligned 0 times
39360173 (64.94%) aligned exactly 1 time
12084173 (19.94%) aligned >1 times
84.88% overall alignment rate
[samopen] SAM header is present: 93 sequences.
SRR1042595.fastq
14603295 reads; of these:
14603295 (100.00%) were unpaired; of these:
918028 (6.29%) aligned 0 times
10403045 (71.24%) aligned exactly 1 time
3282222 (22.48%) aligned >1 times
93.71% overall alignment rate
[samopen] SAM header is present: 93 sequences.
SRR1042596.fastq
65911151 reads; of these:
65911151 (100.00%) were unpaired; of these:
10561790 (16.02%) aligned 0 times
42271498 (64.13%) aligned exactly 1 time
13077863 (19.84%) aligned >1 times
83.98% overall alignment rate
[samopen] SAM header is present: 93 sequences.
SRR1042597.fastq
22210858 reads; of these:
22210858 (100.00%) were unpaired; of these:
1779568 (8.01%) aligned 0 times
15815218 (71.20%) aligned exactly 1 time
4616072 (20.78%) aligned >1 times
91.99% overall alignment rate
[samopen] SAM header is present: 93 sequences.
SRR1042598.fastq
58068816 reads; of these:
58068816 (100.00%) were unpaired; of these:
8433671 (14.52%) aligned 0 times
37527468 (64.63%) aligned exactly 1 time
12107677 (20.85%) aligned >1 times
85.48% overall alignment rate
[samopen] SAM header is present: 93 sequences.
SRR1042599.fastq
24019489 reads; of these:
24019489 (100.00%) were unpaired; of these:
1411095 (5.87%) aligned 0 times
17528479 (72.98%) aligned exactly 1 time
5079915 (21.15%) aligned >1 times
94.13% overall alignment rate
[samopen] SAM header is present: 93 sequences.
SRR1042600.fastq
76361026 reads; of these:
76361026 (100.00%) were unpaired; of these:
8442054 (11.06%) aligned 0 times
50918615 (66.68%) aligned exactly 1 time
17000357 (22.26%) aligned >1 times
88.94% overall alignment rate
[samopen] SAM header is present: 93 sequences.

可以看到比对非常成功!!!我这里就不用表格的形式来展现了,毕竟我又不是给客户写报告,大家就将就着看吧。

05

自学CHIP-seq分析第四讲~必要软件安装以及文章结果下载

博文的顺序有点乱,因为怕读到前面的公共测序数据下载这篇文章的朋友搞不清楚,我如何调用各种软件的,所以我这里强势插入一篇博客来描述这件事,当然也只是略过,我所有的软件理论上都是安装在我的home目录下的biosoft文件夹,所以你看到我一般安装程序都是:

cd ~/biosoft
mkdir macs2 && cd macs2 ##指定的软件安装在指定文件夹里面 Continue reading

05

自学CHIP-seq分析第三讲~公共测序数据下载

这一步跟自学其它高通量测序数据处理一样,就是仔细研读paper,在里面找到作者把原始测序数据放在了哪个公共数据库里面,一般是NCBI的GEO,SRA,本文也不例外,然后解析样本数,找到下载链接规律
## step1 : download raw data
cd ~
mkdir CHIPseq_test && cd CHIPseq_test
mkdir rawData && cd rawData
## batch download the raw data by shell script :
很容易就下载了8个测序文件,每个样本的数据大小,测序量如下
621M Jun 27 14:03 SRR1042593.sra (16.9M reads)
2.2G Jun 27 15:58 SRR1042594.sra (60.6M reads)
541M Jun 27 16:26 SRR1042595.sra (14.6M reads)
2.4G Jun 27 18:24 SRR1042596.sra (65.9M reads)
814M Jun 27 18:59 SRR1042597.sra (22.2M reads)
2.1G Jun 27 20:30 SRR1042598.sra (58.1M reads)
883M Jun 27 21:08 SRR1042599.sra (24.0M reads)
2.8G Jun 28 11:53 SRR1042600.sra (76.4M reads)
 虽然下载的SRA格式数据也是一个很流行的标准,但它只是数据压缩的标准,几乎没有软件能直接跟SRA的格式的测序数据来进行分析,我们需要转成fastq格式,代码如下:
## step2 :  change sra data to fastq files.
## cell line: MCF7 //  Illumina HiSeq 2000 //  50bp // Single ends // phred+33
ls *sra |while read id; do ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump $id;done
rm *sra
解压的详情如下,可以看到SRA格式有6~9倍的压缩了,比zip格式压缩的2~3倍高多了
##  621M --> 3.9G
##  2.2G --> 14G
##  541M --> 3.3G
##  2.4G --> 15G
05

自学CHIP-seq分析第二讲~学习资料的搜集

我只能说,CHIP-seq的确是非常完善的NGS流程了,各种资料层出不穷,大家首先可以看下面几个完整流程的PPT来对CHIP-seq流程有个大致的印象,我对前面提到的文献数据处理的几个要点,就跟下面这个图片类似:
chip-seq-workflow-all-5-steps
QuEST is a statistical software for analysis of ChIP-Seq data with data and analysis results visualization through UCSC Genome Browser.  http://www-hsc.usc.edu/~valouev/QuEST/QuEST.html
MeDIP-seq and histone modification ChIP-seq analysis  http://crazyhottommy.blogspot.com/2014/01/medip-seq-and-histone-modification-chip.html
一个实际的CHIP-seq数据分析例子: http://www.biologie.ens.fr/~mthomas/other/chip-seq-training/

Continue reading

05

自学CHIP-seq分析第一讲~文献选择与解读

文章:CARM1 Methylates Chromatin Remodeling Factor BAF155 to Enhance Tumor Progression and Metastasis 

我很早以前想自学CHIP-seq的时候就关注过这篇文章,那时候懂得还不多,甚至都没有仔细看这篇文章就随便下载了数据进行分析,也只是跑一些软件而已,这次仔细阅读这篇文章才发现里面的门道很多,尤其是CHIP-seq的实验基础,以及表观遗传学的生物学基础知识,我有时间一定要把这篇文章翻译一下。学习这篇文章前一定要温习一些生物学知识,见我上一篇博客

Continue reading

05

生物学基础知识~CARM和SWI/SNF复合物

因为最近在研究CHIP-seq测序数据处理,发现有些文章重点并不是数据处理本身,而是对生物学基础知识的掌控以及实验设计,这篇文章我重点推荐一下,我略微翻译了一些,笔记如下:

SWI/SNF(BAF) chromatin remodeling complex  染色质重构复合物,以及被广泛发现在各种癌症患者体内均有突变,这个复合物利用ATP水解释放的能量来驱动核小体运动以及调控染色质的结构。
发现历史: 最新是在yeast里面发现它的突变会影响 mating type SWItching 和 导致sucrose nonfermenting 的表型,所以才简称为SWI/SNF,在哺乳动物体内也被广泛研究现在,它被发现参与了细胞分化、增殖、还有各种DNA修复功能的实现。
组成结构:本质是一个蛋白质复合物,约15个亚基,每个亚基都由一个独立的基因转录翻译而来,每个基因都有专门的文章研究过。
SWI-SNF-geneFamily-SMARCC1-BAF155
生物功能:早在1998年就被发现与癌症相关,在GO数据库里面还可以查到相关资料,所以很容易在各种通路分析结果里面看到它的身影
GO-0016514-SWI-SNF-complex
随着癌症基因组测序的进展,至少有8种染色质重构复合物的亚基有recurrent mutation情况,然后有一句话说得特别好:
This has resulted in interest in indentifying mechanisms by which activity of the SWI/SNF complex is regulated, with the hope that such mechanistic understanding may reveal novel opportunities for therapeutic intervention .
CARM1基因是PRMT系列的一种:protein arginine Methyltransferases一直与基因的转录与翻译相关,主要功能就是使得蛋白质的arginine甲基化,这个基因家族目前有9个基因,只有CARM1命名比较奇葩。
PRMT-gene-family-CARM1
其中RM都是 Arginine Methyltransferase的简称,很容易理解,
CARM1这个酶的底物不仅仅包括参与染色质重构,还有基因转录调控,剪切因子,组蛋白乙酰化还有RNA结合蛋白
chromatin remodeling (  SWI/SNF(BAF) chromatin remodeling complex  )
gene transcription   (histone H3(at R17))
histone acetyltransferases  ( p300/CBP )
splicing factors  ( CA150,SAP49,SmB,U1C  )
RNA-binding proteins   ( PABP1, HuR, HuD)
也有实验证明非常多的癌症种类病人都发现了CARM1表达异常升高,而且刺激乳腺癌的恶化,而且是很多癌症相关恶化蛋白的共刺激因子,比如p53,E2F1,NF-kB,b-catenin,steroid hormone receptors. 
但是CARM1到底是如何在乳腺癌体内发生作用的,其中机制并不完全清楚。
04

跟师妹聊Exome-seq、ChIP-seq、RNA-seq之间的差异

最近学习CHIP-seq的分析流程,略有点心得,也跟以前掌握的WES和RNA-seq做了一些比较,趁跑步的时候跟师妹讨论了一下,正好师妹写了一篇博客来分享这个讨论结果,我也借此机会转载过来,分享给大家,算是借花献佛吧!师妹的博文是用markdown写作,我觉得大家应该直接看她的文章,写得条理清楚:Exome-seq、ChIP-seq、RNA-seq之间的差异 Continue reading

03

自学miRNA-seq分析第八讲~miRNA-mRNA表达相关下游分析

通过前面的分析,我们已经量化了ET1刺激前后的细胞的miRNA和mRNA表达水平,也通过成熟的统计学分析分别得到了差异miRNA和mRNA,这时候我们就需要换一个参考文献了,因为前面提到的那篇文章分析的不够细致,我这里选择了浙江大学的一篇TCGA数据挖掘分析文章Identifying miRNA/mRNA negative regulation pairs in colorectal cancer,里面首先就是查找miRNA-mRNA基因对,因为miRNA主要还是负向调控mRNA表达,所以根据我们得到的两个表达矩阵做相关性分析,很容易得到符合统计学意义的miRNA-mRNA基因对,具体分析内容如下:

把得到的差异miRNA的表达量画一个热图,看看它是否能显著的分类
用miRWalk2.0等数据库或者根据来获取这些差异miRNA的validated target genes
然后看看这些pairs of miRNA- target genes的表达量相关系数,选取显著正相关或者负相关的pairs
这些被选取的pairs of miRNA- target genes拿去做富集分析
最后这些pairs of miRNA- target genes做PPI网络分析

首先我们看第一个热图的实现:

resOrdered=na.omit(resOrdered)
DEmiRNA=resOrdered[abs(resOrdered$log2FoldChange)>log2(1.5) & resOrdered$padj <0.01 ,]
write.csv(resOrdered,"deseq2.results.csv",quote = F)
DEmiRNAexprSet=exprSet[rownames(DEmiRNA),]
write.csv(DEmiRNAexprSet,'DEmiRNAexprSet.csv')

DEmiRNAexprSet=read.csv('DEmiRNAexprSet.csv',stringsAsFactors = F)
exprSet=as.matrix(DEmiRNAexprSet[,2:7])
rownames(exprSet)=rownames(DEmiRNAexprSet)
heatmap(exprSet)
gplots::heatmap.2(exprSet)
library(pheatmap)
## http://biit.cs.ut.ee/clustvis/

因为我前面保存的表达量就基于counts的,所以画热图还需要进行normalization,我这里懒得弄了,就用了一个网页版工具,自动出热图http://biit.cs.ut.ee/clustvis/

miRNA-heatmap

感觉还不错,可以很清楚的看到ET1刺激前后细胞中miRNA表达量变化

然后就是检验我们选取的感兴趣的有显著差异的miRNA的target genes,这时候有两种方法,一个是先由数据库得到已经被检验的miRNA的target genes,另一种是根据miRNA和mRNA表达量的相关性来预测。

用数据库来查找MiRNA的作用基因,非常多的工具,比较常用的有TargetScan/miRTarBase
### http://nar.oxfordjournals.org/content/early/2015/11/19/nar.gkv1258.full
### http://mirtarbase.mbc.nctu.edu.tw/
### http://mirtarbase.mbc.nctu.edu.tw/cache/download/6.1/hsa_MTI.xlsx
### http://www.targetscan.org/vert_71/ (version 7.1 (June 2016))
我还看到过一个整合工具: miRecords  (DIANA-microT, MicroInspector, miRanda, MirTarget2, miTarget, NBmiRTar, PicTar, PITA, RNA22, RNAhybrid and TargetScan/TargertScanS)里面提到了查找MiRNA的作用基因这一过程,高假阳性,至少被5种工具支持,才算是真的
还有很多类似的工具,miRWalk2,psRNATarget网页版工具,最后值得一提的是中山大学的: starBase  Pan-Cancer Analysis Platform is designed for deciphering Pan-Cancer Networks of lncRNAs, miRNAs, ceRNAs and RNA-binding proteins (RBPs) by mining clinical and expression profiles of 14 cancer types (>6000 samples) from The Cancer Genome Atlas (TCGA) Data Portal (all data available without limitations).虽然我没有仔细的用,但是看介绍好牛的样子,还有一个R包:miRLAB我玩了一会,它是先通过算所有配对的miRNA- genes的表达量相关系数,选取显著正相关或者负相关的pairs,然后反过来通过已知数据库来验证。

后面我就不讲了,主要看你得到miRNA的时候其它生物学数据是否充分,如果是癌症病人,有生存相关数据,可以做生存分析,如果你同时测了甲基化数据,可以做甲基化相关分析~~~~~~~~~

如果只是单纯的miRNA测序数据,可以回过头去研究一下de novo的miRNA预测的步骤,也是研究重点

 

01

自学miRNA-seq分析第七讲~miRNA样本配对mRNA表达量获取

这一讲其实算不上是自学miRNA-seq分析,本质就是affymetrix的mRNA表达芯片数据分析,而且还是最常用的那种GPL570    HG-U133_Plus_2,但是因为是跟miRNA样本配对检测的,而且后面会利用到这两个数据分析结果来做共表达网络分析等等,所以就贴出对该芯片数据的分析结果。文章里面也提到了 Messenger RNA expression analysis identified 731 probe sets with significant differential expression,作者挑选的差异分析结果的显著基因列表如下: Continue reading

01

自学miRNA-seq分析第六讲~miRNA表达量差异分析

这一讲是miRNA-seq数据分析的分水岭,前面的5讲说的是读文献下载数据比对然后计算表达量,属于常规的流程分析,一般在公司测序之后都可以拿到分析结果,或者文献也会给出下载结果。但是单纯的分析一个样本意义不大,一般来说,我们做研究都是针对于不同状态下的miRNA表达量差异分析,然后做注释,功能分析,网络分析,这才是重点,也是难点。我这里就直接拿文献处理好的miRNA表达量来展示如何做下游分析,首先就是差异分析啦: Continue reading