30

找变异的流程

找变异简单点说,就是把高通量测序得到的成千上万条序列片段比对到合适的参考基因组,找到那些成

功比对的片段与参考基因组的微小差异情况。 那么就涉及到存储测序数据的fastq数据格式,比对的工具,比对后的sam格式,找微小差异的工具,差异结果的vcf文件,每个步骤的软件选择,参数 调整。当然,最重要的是走通整个流程,明白自己在做什么。

Continue reading

21

自学无参RNAseq数据分析第一讲之参考文献解读

这是我为新创办的 生信技能树 论坛写的帖子,也适合本博客,所以转载过来: http://www.biotrainee.com/thread-243-1-1.html 

以前做的都是有参转录组分析,只需要找到参考基因组和注释文件,然后走QC-->alignment-->counts->DEG-->annotation的流程图即可。
现在开始学习新的东西了,就是无参转录组分析,这里记录一下自己的学习笔记,首先还是资料收集,这次,我就针对性的看5个 全流程化的转录组 de novo 分析 文章,如下:
http://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-15-554  2014年栀子花的花瓣衰老的标准de novo 转录组分析,数据如下:用Trinity做组装,用NCBI non-redundant (Nr) database库做注释,做了差异分析(栀子花花期分成4个阶段),GO/KEGG注释,然后做了RT-qPCR的实验验证。
多做了一个 Clusters of Orthologus Groups (COG)的数据库注释

Raw Reads
Clean Reads
Contigs
Unigenes
Annotated
Transcriptome
55,092,396
50,335,672
102,263
57,503
39,459

 

http://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-15-236  2014 巴西橡胶树的研究,是一个综合多组织样本的RNA库,ployT建库,454测序,用的是est2Assembly 和gsassembler 软件做组装,用 NCBI RefSeq, Plant Protein Database 做注释,因为没有分组,所以不必做差异分析,只需要找SNV和SSR标记即可,最后也是做GO/KEGG注释

https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-016-2633-2 2015 萝卜,用illumina进行转录组测序,用Trinity组装,用RPKM值算unigene的表达量,也是用 BLASTx来对Trinity结果进行注释,注释到NR,NT,Swiss-Prot,GO,COG,kegg数据库,其中GO注释用的是Blast2GO,最后也做了RT-qPCR 实验验证,某些基因在leaf里面的表达量显著高于其它tissue,有原始数据:http://www.ncbi.nlm.nih.gov/sra/?term=SRX1671013
转录组分析结果结果:A total of 54.64 million clean reads and 111,167 contigs representing 53,642 unigenes were obtained from the radish leaf transcriptome.

http://www.nature.com/articles/srep08259 2015 芹菜 叶片发育中木质素的探究,测序的reads是A total of 32,477,416 quality reads were recorded for the leaves at Stage 1, 53,675,555 at Stage 2, and 27,158,566 at Stage 3, respectively.,也是用Trinity组装,kmer值设为25,组装结果:33,213 unigenes with an average length of 1,478 bp, a maximum length of 17,075 bp, and an N50 of 2,060 bp,然后用eggNOG/GO/KEGG数据库来注释。文章正文给了所用到的软件和数据库的详细链接
最后还用了 real-time PCR assays          来看 roots, stems, petioles, and leaf blade 这些组织的基因表达差异情况

http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0128659 对 三疣梭子蟹 的卵巢和睾丸的转录组研究,,也是标准的转录组de novo 分析流程,非常值得借鉴
NCBI有上传原始数据:SRR1920180  和SRR1920180  

总结好这5篇文献的数据分析流程,就差不多明白如何做无参的转录组de novo分析了

09

RNAseq数据完整生物信息分析流程第一讲之文献数据下载

我这里拿的是bioconductor里面最常用的airway数据,因为差异表达分析在bioconductor里面是重点,它们这些包在介绍自己的算法以及做示范的时候都用的这个数据。可以在GEO数据库里面看到信息描述:http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52778  可以看到是Illumina HiSeq 2000 (Homo sapiens) ,75bp paired-end 这个信息很重要,决定了下载sra数据之后如何解压以及如何比对。也可以看到作者把所有的测序原始数据都上传到了SRA中心:http://www.ncbi.nlm.nih.gov/sra?term=SRP033351  ,这里可以在linux服务器上面写一个简单的脚本批量下载所有的测序数据,然后根据GEO里面描述的metadata把原始数据改名。

for ((i=508;i<=523;i++)) ;do wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP033/SRP033351/SRR1039$i/SRR1039$i.sra;done
ls *sra |while read id; do ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 $id;done

需要自己看SRA里面的数据记录,上面的脚本不难写出,然后因为是Illumina的双端测序,所以我们用fastq-dump --split-3命令来把sra格式数据转换为fastq,但是因为这里有16个测序数据,所以最好是同步改名,我这里用脚本批量生成改名脚本如下:

为了节省空间,我用了--gzip压缩,该文件名,用-A参数。

nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N61311_untreated SRR1039508.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N61311_Dex SRR1039509.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N61311_Alb SRR1039510.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N61311_Alb_Dex SRR1039511.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N052611_untreated SRR1039512.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N052611_Dex SRR1039513.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N052611_Alb SRR1039514.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N052611_Alb_Dex SRR1039515.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N080611_untreated SRR1039516.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N080611_Dex SRR1039517.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N080611_Alb SRR1039518.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N080611_Alb_Dex SRR1039519.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N061011_untreated SRR1039520.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N061011_Dex SRR1039521.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N061011_Alb SRR1039522.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N061011_Alb_Dex SRR1039523.sra &

可以看到这里的16个样本来源于同样的4个人,是HASM细胞系,处理详情如下:

测序基础:
HASM细胞系-human airway smooth muscle,
The Illumina TruSeq assay was used to prepare 75bp paired-end libraries for HASM cells from four white male donors under four treatment conditions:
1) no treatment;
2) treatment with a β2-agonist (i.e. Albuterol, 1μM for 18h);
3) treatment with a glucocorticosteroid (i.e. Dexamethasone (Dex), 1μM for 18h);
4) simultaneous treatment with a β2-agonist and glucocorticoid
and the libraries were sequenced with an Illumina Hi-Seq 2000 instrument.
我们这里只是先根据fastq数据比对到参考基因组,然后计算每个样本的表达量即可,后续的分组计算差异表达,就需要个性化了。

下载的sra大小如下:

-rw-rw-r-- 1 jmzeng jmzeng 1.6G Aug 9 04:21 SRR1039508.sra
-rw-rw-r-- 1 jmzeng jmzeng 1.5G Aug 9 05:20 SRR1039509.sra
-rw-rw-r-- 1 jmzeng jmzeng 1.6G Aug 9 06:14 SRR1039510.sra
-rw-rw-r-- 1 jmzeng jmzeng 1.5G Aug 9 07:05 SRR1039511.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.1G Aug 9 08:07 SRR1039512.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.3G Aug 9 09:17 SRR1039513.sra
-rw-rw-r-- 1 jmzeng jmzeng 3.1G Aug 9 10:56 SRR1039514.sra
-rw-rw-r-- 1 jmzeng jmzeng 1.9G Aug 9 11:56 SRR1039515.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.1G Aug 9 13:02 SRR1039516.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.6G Aug 9 14:16 SRR1039517.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.3G Aug 9 15:17 SRR1039518.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.0G Aug 9 16:05 SRR1039519.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.1G Aug 9 16:56 SRR1039520.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.4G Aug 9 17:57 SRR1039521.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.0G Aug 9 18:46 SRR1039522.sra
-rw-rw-r-- 1 jmzeng jmzeng 1.4G Aug 9 19:28 SRR1039523.sra

解压后成双端测序的fastq数据如下:

 -rw-rw-r-- 1 jmzeng jmzeng 2.5G Aug 9 20:12 N052611_Alb_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 2.5G Aug 9 20:12 N052611_Alb_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 20:44 N052611_Alb_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 20:44 N052611_Alb_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 289M Aug 9 20:44 N052611_Alb_Dex.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 951M Aug 9 20:59 N052611_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 954M Aug 9 20:59 N052611_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.7G Aug 9 20:53 N052611_untreated_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.7G Aug 9 20:53 N052611_untreated_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.5G Aug 9 20:45 N061011_Alb_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.5G Aug 9 20:45 N061011_Alb_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.9G Aug 9 20:59 N061011_Alb_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.9G Aug 9 20:59 N061011_Alb_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 16M Aug 9 20:45 N061011_Alb.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.4G Aug 9 20:48 N061011_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.4G Aug 9 20:48 N061011_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.2G Aug 9 20:00 N061011_untreated_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.2G Aug 9 20:00 N061011_untreated_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 759M Aug 9 20:00 N061011_untreated.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.9G Aug 9 20:03 N080611_Alb_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.9G Aug 9 20:03 N080611_Alb_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 19:59 N080611_Alb_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 19:59 N080611_Alb_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 535M Aug 9 19:59 N080611_Alb_Dex.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 2.1G Aug 9 20:06 N080611_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 2.1G Aug 9 20:06 N080611_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.6G Aug 9 20:01 N080611_untreated_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.6G Aug 9 20:01 N080611_untreated_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:09 N61311_Alb_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:09 N61311_Alb_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:08 N61311_Alb_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:08 N61311_Alb_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.2G Aug 9 08:07 N61311_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.2G Aug 9 08:07 N61311_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:09 N61311_untreated_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:09 N61311_untreated_2.fastq.gz

接下来所有的分析就基于此数据啦

 

 

02

生信人必学ftp站点之1000genomes

千人基因组计划的重要性我也不想多说了,由于时间跨度比较长,最终的数据不只是一千人,最新版共有NA编号开头的1182个人,HG开头的1768个人!它的官方网站是:有一个ppt讲得很清楚如何通过官网做的data portal来下载数据:https://www.genome.gov/pages/research/der/ichg-1000genomestutorial/how_to_access_the_data.pdf 我不喜欢可视化的界面,我比较喜欢直接进入ftp自己翻需要的数据,千人基因组计划不仅仅有自己的ftp站点,而且在NCBI,EBI和sanger研究所里面也有数据源可以下载, 是非常丰富的生信入门资源!

Continue reading

02

生信人必学ftp站点之NCBI-GEO

NCBI的重要性我就不多说了,Gene Expression Omnibus database (GEO)是由NCBI负责维护的一个数据库,设计初衷是为了收集整理各种表达芯片数据,但是后来也加入了甲基化芯片,lncRNA,miRNA,CNV芯片等各种芯片,甚至高通量测序数据!所有的数据均可以在ftp站点下载:ftp://ftp-trace.ncbi.nih.gov/geo/ Continue reading

14

ChIP-Seq文献数据重新分析解读第二例

这篇文章是朋友推荐的, 我觉得作为CHIP-seq学习材料再好不过了,所以推荐给大家。是全基因组范围的BRCA1和PALB2的转录共激活机制的探究。请务必先看我的CHIP-seq自学系列教程,跟着好好学习!数据如下:
GSM997540    BRCA1    SRR553473.sra    Read 18878514 spots
GSM997541    PALB2    SRR553474.sra    Read 17615498 spots
GSM997542    P_Ser2    SRR553475.sra    Read 35396009 spots

Continue reading

13

ChIP-Seq文献数据重新分析解读第一例

文章是:Genome-wide maps of H3K4me2/3 in prostate cancer cell line LNCaP,数据在GEO可以下载。GSE20042,下面的所有分析,需要26G的空间。
作者想看看用 dihydrotestosterone (雄激素)处理了 cancer cell line LNCaP 这个细胞系之后,看看组蛋白甲基化修饰变化,主要是看H3K4me2和H3K4me3这两种组蛋白甲基化区别,分成三组,分别是处理前,处理后4H和16H,共有5个条件的数据,但是有7个fastq文件。
08

自学lncRNA-seq数据分析第一讲~学习大纲

lncRNA分析跟常见的mRNA-seq分析重合度很高,无非也是把测序的fastq文件mapping到参加基因组,获取转录本信息,转录本表达定量,表达量的差异分析,比较新的分析就是把转录本分成了lncRNA和mRNA,这样可以考虑它们之间的互相作用,也可以在实验设计的时候加入miRNA和CHIP-seq,这样多种数据结合分析,显得更高大上一点,也能更好的刻画机体状态,从而回答生物学假设,我这里先列出我的自学大纲,如果有朋友想跟着我学习,可以发email给我,我的邮箱是jmzeng1314,是163邮箱,发信给我前,先看如果你希望我回答你的问题 这篇文章:

Continue reading

07

自学CHIP-seq分析第九讲~CHIP-seq可视化大全

讲到这里,我们的自学CHIP-seq分析系列教程就告一段落了,当然,我会随时查漏补缺,根据读者的反馈来更新着系列教程。其实可视化这已经是一个比较复杂的方向了,不仅仅是针对于CHIP-seq数据。可视化本身是发文章的先决条件,而让人一目了然图片也说明了数据分析人员对数据本身的理解。我这里就列出一些目录和一些工具,和ppt。这个主要靠大家自学了,而且我博客空间有限,就不上传一大堆图片了,大家随便找一些经典的paper里面都会有很多可视化分析。

首先强烈推荐两个网页版工具,针对找到的peaks可视化:

然后再推荐一个哈佛刘小乐实验室出品的软件,也是专门为了作图http://liulab.dfci.harvard.edu/CEAS/usermanual.html

Continue reading

07

自学CHIP-seq分析第八讲~寻找motif

motif是比较有特征的短序列,会多次出现的,一般认为它的生物学意义重大,做完CHIP-seq分析之后,一般都会寻找motif 。查找有两种,一种是de novo的,要求的输入文件的fasta序列,一般是根据peak的区域的坐标提取好序列 。另一种是依赖于数据库的搜寻匹配,很多课题组会将现有的ChIP-seq数据进行整合,提供更全面,更准确的motif数据库。

Continue reading

06

自学CHIP-seq分析第七讲~peaks注释

经过前面的CHIP-seq测序数据处理的常规分析,我们已经成功的把测序仪下机数据变成了BED格式的peaks记录文件,我选取的这篇文章里面做了4次CHIP-seq实验,分别是两个重复的野生型MCF7细胞系的 BAF155 immunoprecipitates和两个重复的突变型MCF7细胞系的 BAF155 immunoprecipitates,这样通过比较野生型和突变型MCF7细胞系的 BAF155 immunoprecipitates的结果的不同就知道该细胞系的BAF155 突变,对它在全基因组的结合功能的影响啦。

#我这里直接从GEO里面下载了peaks结果,它们详情如下:wc -l *bed
6768 GSM1278641_Xu_MUT_rep1_BAF155_MUT.peaks.bed
3660 GSM1278643_Xu_MUT_rep2_BAF155_MUT.peaks.bed
11022 GSM1278645_Xu_WT_rep1_BAF155.peaks.bed
5260 GSM1278647_Xu_WT_rep2_BAF155.peaks.bed
49458 GSM601398_Ini1HeLa-peaks.bed
24477 GSM601398_Ini1HeLa-peaks-stringent.bed
12725 GSM601399_Brg1HeLa-peaks.bed
12316 GSM601399_Brg1HeLa-peaks-stringent.bed
46412 GSM601400_BAF155HeLa-peaks.bed
37920 GSM601400_BAF155HeLa-peaks-stringent.bed
30136 GSM601401_BAF170HeLa-peaks.bed
25432 GSM601401_BAF170HeLa-peaks-stringent.bed

每个BED的peaks记录,本质是就3列是需要我们注意的,就是染色体,以及在该染色体上面的起始和终止坐标,如下:

#PeakID chr start end strand Normalized Tag Count region size findPeaks Score Clonal Fold Change
chr20 52221388 52856380 chr20-8088 41141 +
chr20 45796362 46384917 chr20-5152 31612 +
chr17 59287502 59741943 chr17-2332 29994 +
chr17 59755459 59989069 chr17-667 19943 +
chr20 52993293 53369574 chr20-7059 12642 +
chr1 121482722 121485861 chr1-995 9070 +
chr20 55675229 55855175 chr20-6524 7592 +
chr3 64531319 64762040 chr3-4022 7213 +
chr20 49286444 49384563 chr20-4482 6165 +

我们所谓的peaks注释,就是想看看该peaks在基因组的哪一个区段,看看它们在各种基因组区域(基因上下游,5,3端UTR,启动子,内含子,外显子,基因间区域,microRNA区域)分布情况,但是一般的peaks都有近万个,所以需要批量注释,如果脚本学的好,自己下载参考基因组的GFF注释文件,完全可以自己写一个,我这里会介绍一个R的bioconductor包ChIPpeakAnno来做CHIP-seq的peaks注释,下面的包自带的示例:

library(ChIPpeakAnno)
bed <- system.file("extdata", "MACS_output.bed", package="ChIPpeakAnno")
gr1 <- toGRanges(bed, format="BED", header=FALSE)
## one can also try import from rtracklayer
library(rtracklayer)
gr1.import <- import(bed, format="BED")
identical(start(gr1), start(gr1.import))
gr1[1:2]
gr1.import[1:2] #note the name slot is different from gr1
gff <- system.file("extdata", "GFF_peaks.gff", package="ChIPpeakAnno")
gr2 <- toGRanges(gff, format="GFF", header=FALSE, skip=3)
ol <- findOverlapsOfPeaks(gr1, gr2)
makeVennDiagram(ol)

##还可以用binOverFeature来根据特定的GRanges对象(通常是TSS)来画分布图
## Distribution of aggregated peak scores or peak numbers around transcript start sites.

可以看到这个包使用起来非常简单,只需要把我们做好的peaks文件(GSM1278641_Xu_MUT_rep1_BAF155_MUT.peaks.bed等等)用toGRanges或者import读进去,成一个GRanges对象即可,上面的代码是比较两个peaks文件的overlap。然后还可以根据R很多包都自带的数据来注释基因组特征:

data(TSS.human.GRCh37) ## 主要是借助于这个GRanges对象来做注释,也可以用getAnnotation来获取其它GRanges对象来做注释
## featureType : TSS, miRNA, Exon, 5'UTR, 3'UTR, transcript or Exon plus UTR
peaks=MUT_rep1_peaks
macs.anno <- annotatePeakInBatch(peaks, AnnotationData=TSS.human.GRCh37,
output="overlapping", maxgap=5000L)

## 得到的macs.anno对象就是已经注释好了的,每个peaks是否在基因上,或者距离基因多远,都是写的清清楚楚
if(require(TxDb.Hsapiens.UCSC.hg19.knownGene)){
aCR<-assignChromosomeRegion(peaks, nucleotideLevel=FALSE,
precedence=c("Promoters", "immediateDownstream",
"fiveUTRs", "threeUTRs",
"Exons", "Introns"),
TxDb=TxDb.Hsapiens.UCSC.hg19.knownGene)
barplot(aCR$percentage)
}

得到的条形图如下,虽然很丑,但这就是peaks注释的精髓,搞清楚每个peaks在基因组的位置特征:

ChIPpeakAnno-genomic-region-distribution

 

同理,对每个peaks文件,都可以做类似的分析!

但是对多个peaks文件,比如本文中的,想比较野生型和突变型MCF7细胞系的 BAF155 immunoprecipitates的结果的不同,就需要做peaks之间的差异分析,已经后续的差异基因注释啦

当然,值得一提的是peaks注释我更喜欢网页版工具,反正peaks文件非常小,直接上传到别人做好的web tools,就可立即出一大堆可视化图表分析结果啦,大家可以去试试看:

虽然我花费了大部分篇幅来描述ChIPpeakAnno这个包的用法,但是真正的重点是你得明白peaks记录了什么,要注释什么,已经把这3个网页工具的可视化图表分析结果全部看懂,这网页版工具才是重点!!!
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

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