06

Bioconductor包chimeraviz嵌合RNA可视化

Bioconductor包chimeraviz嵌合RNA可视化

高通量RNA测序已经能够更高效地检测融合转录本,但是融合检测的技术和相关软件通常产生高错误发现率。而一个自动整合RNA数据和已知基因组特征的可视化框架对于结果的检验是有帮助的。2017年发布的一个bioconductor包,chimeraviz就可以做到自动创建嵌合RNA可视化。

支持来自9种不同融合发现工具(deFuseEricScript、InFusion、JAFFA、FusionCatcher、FusionMap、PRADA、SOAPfuse和STAR-FUSION)的输入。 Continue reading

05

用LeafCutter探索转录组数据的可变剪切

用LeafCutter探索转录组数据的可变剪切

该软件早在2016年就公布了,发表在biorxiv预印本上面,但直到2017年的双11,才发表在NG上面,文章是 : Annotation-free quantification of RNA splicing using LeafCutter 最大的特点应该是不需要参考基因组的基因注释信息了吧,就是gtf/gff文件可以省略,当然,比对还是需要的。它还有另外一个非常重要的功能,splicing quantitative trait loci (sQTLs) 但是跟我目前关系不大, 就不介绍了。 Continue reading

十二 14

使用SGSeq探索可变剪切

可变剪切是指mRNA前体以多种方式将exon连接在一起的过程。 由于可变剪切使一个基因产生多个mRNA转录本,不同mRNA可能翻译成不同蛋白。

可变剪切背景知识

转录组一般是指从细胞或组织的基因组所转录出来的RNA的总和,包括编码蛋白质的mRNA和各种非编码RNA(rRNA,tRNA,snRNA,snoRNA,lncRNA,microRNA等)。真核生物的基因结构是不连续的,如下图:

Continue reading

十一 02

一个植物转录组项目的实战

其实这个植物是拟南芥,所以跟人类研究的数据处理大同小异。

转录组

转录组测序的研究对象为特定细胞在某一功能状态下所能转录出来的所有 RNA 的总和,包括 mRNA 和非编码 RNA 。通过转录组测序,能够全面获得物种特定组织或器官的转录本信息,从而进行转录本结构研究、变异研究、基因表达水平研究以及全新转录本发现等研究。 Continue reading

十二 30

一个RNA-seq实战-超级简单-2小时搞定!

请不要直接拷贝我的代码,需要自己理解,然后打出来,思考我为什么这样写代码。
软件请用最新版,尤其是samtools等被我存储在系统环境变量的,考虑到读者众多,一般的软件我都会自带版本信息的!
我用两个小时,不代表你是两个小时就学会,有些朋友反映学了两个星期才 学会,这很正常,没毛病,不要异想天开两个小时就达到我的水平。
转录组如果只看表达量真的是超级简单,真是超级简单,而且人家作者本来就测是SE50,这种破数据,也就是看表达量用的!
首先作者分析结果是:
1

Continue reading

十一 25

hisat2+stringtie+ballgown

早在去年九月,我就写个博文说 RNA-seq流程需要进化啦! http://www.bio-info-trainee.com/1022.html  ,主要就是进化成hisat2+stringtie+ballgown的流程,但是我一直没有系统性的讲这个流程,因为我觉真心木有用。我只用了里面的hisat来做比对而已!但是群里的小伙伴问得特别多,我还是勉为其难的写一个教程吧,你们之间拷贝我的代码就可以安装这些软件的!然后自己找一个测试数据,我的脚本很容易用的! Continue reading

十一 15

htseq-counts跟bedtools的区别

我以前写过bedtools和htseq-counts的教程,它们都可以用来对比对好的bam文件进行计数,正好群里有小伙伴问我它们的区别,我就简单做了一个比较,大家可以先看看我以前写的软件教程。写的有的挫:

使用Bedtools对RNA-seq进行基因计数 ,

转录组HTseq对基因表达量进行计数

言归正传,我这里没精力去探究它们的具体原理,只是看看它们数一个read是否属于某个基因的时候,区别在哪里,大家看下图: Continue reading

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

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

 

 

12

转录组 de novo流程–包括转录本完整注释

有网友咨询过对于没有参考基因组或者转录组的物种,如何做RNA-seq分析。我觉得这个问题太大了,而且我还真的对这个没有经验。但是我以前看到过一篇文献,里面提到过一个非常全面的转录组 de novo组装注释流程,所以我摘抄了文章里面的生物信息学处理部分,分享给大家: Continue reading

15

融合基因检测软件-soapfusion

开发单位:华大,SOAP系列软件套装!

功能:检测合基因
优点:在现有的各种软件里面表现算是最好的
算法:是hash index,跟其它bwt算法不太一样
其它软件有: FusionSeq [21], deFuse [22], TopHat-Fusion [23], FusionHunter [24], SnowShoes-FTD [25], chimerascan [26] and FusionMap [27]
具体的算法我没看,因为只是有需求,正好有一些RNA-seq数据又想看看样本融合基因情况。所以就测试这个软件,通俗点说,融合基因原理其实很简单,如果有足够多的reads一部分比对到一个基因,另一部分比对到另一个基因,就可以说明它们两个基因发生了融合现象!如果是PE测序,那么更方便,左右两端reads比对情况也可以考虑。我就不多说废话了,直接上教程吧!
一,软件安装
下载压缩包,解压后即可使用!!!
推荐用最新版,然后看作者说明书的时候也要看清楚!
我反正好几次都搞糊涂了,最后联系了作者才搞明白,作者说他想更新到2.0版本,直接用HISAT的比对sam文件来做,但是还在筹备中,我觉得有点悬!
1
解压后是一堆perl程序,都在source目录下,source目录下面还有bin下面附带了几个第三方软件,包括bwa,blast和soap,最后都用得着!
有个很重要的问题,一定要软件自带的perl模块添加到perl的环境变量。不然那些perl程序运行会报错!
配置文件需要修改,就把几个目录放进去即可

二,输入数据准备

这里最重要的就是制作数据库!!!
作者给了非常详细的制作过程,我觉得还是不够清楚,所以再讲一遍!
首先下载5个文件:
6.5K Jun 15  2009 cytoBand.txt.gz
3.0G Oct 12  2012 hg19.fa
2.5M Mar 15 10:30 HGNC_Gene_Family_dataset
38M Feb  8  2014 Homo_sapiens.GRCh37.75.gtf.gz
202 Jan 19 16:07 HumanRef_refseg_symbols_relationship.list

文件下载地址,作者已经给出了!

我把这些文件都放在的当前文件夹下面的raw这个子文件夹,因为我要当前文件夹作为该软件的database文件夹!!!
然后运行命令!
我在SOAPfuse-v1.27文件下面运行:
perl ../SOAPfuse-v1.27/source/SOAPfuse-S00-Generate_SOAPfuse_database.pl  \
-wg raw/hg19.fa  -gtf raw/Homo_sapiens.GRCh37.75.gtf.gz  -cbd raw/cytoBand.txt.gz   -gf raw/HGNC_Gene_Family_dataset \
-rft raw/HumanRef_refseg_symbols_relationship.list \
 -sd ../SOAPfuse-v1.27 -dd ./

这一步耗时很长,4~6小时,创造了transcript.fa和gene.fa,然后还对他们建立bwa和soap的index,所以有点慢!

构建成功会有提示:
Congratulations!
You have constructed SOAPfuse database files successfully.
These database files are all stored in directory you supplied:
/home/jmzeng/biosoft/SOAPfuse/db_v1.27/
They are all generated based on public data files you supplied:
whole_genome_fasta_file:   /home/jmzeng/biosoft/SOAPfuse/db_v1.27/raw/hg19.fa
gtf_annotation_file:       /home/jmzeng/biosoft/SOAPfuse/db_v1.27/raw/Homo_sapiens.GRCh37.75.gtf.gz
Chr_Bandregion_file:       /home/jmzeng/biosoft/SOAPfuse/db_v1.27/raw/cytoBand.txt.gz
HGNC_gene_family_file:     /home/jmzeng/biosoft/SOAPfuse/db_v1.27/raw/HGNC_Gene_Family_dataset
gtf_segname2refseg_list:   /home/jmzeng/biosoft/SOAPfuse/db_v1.27/raw/HumanRef_refseg_symbols_relationship.list
这些目录很重要,接下来制作配置文件会用得着!
To use these database files, just set the 'DB_db_dir' in config file as belowed:
DB_db_dir  =   /home/jmzeng/biosoft/SOAPfuse/db_v1.27
配置文件需要修改下面5个
DB_db_dir = /DATABASE_DIR/
PG_pg_dir = /TOOL_DIR/source/bin
PS_ps_dir = /TOOL_DIR/source
PD_all_out = /out_directory/
PA_all_fq_postfix = PostFix
其实你仔细阅读了说明书,你就知道该修改成什么样子了!
最后制作sample list文件
我这里只有一个sample,所以文件就一句话即可
test test test 100
所以我的有下面两个文件,都是为了顺应作者的需求我才搞了test/test/test这么无聊的东西!!!
/home/jmzeng/test_for_soapfuse/test/test/test_1.fq.gz
/home/jmzeng/test_for_soapfuse/test/test/test_2.fq.gz
如果你有多个sample需要一起运行,你就要仔细读作者的readme了,它把这个配置文件搞得特别复杂!!!

三,运行命令

如果文件都准备好了,运行命令非常简单!!
perl SOAPfuse-RUN.pl -c <config_file> -fd <WHOLE_SEQ-DATA_DIR> -l <sample_list> -o <out_directory> [Options]

运行的非常慢!!!

因为需要重新比对,知道

四,数据结果解读

结果,作者已经说的很清楚了,我就不多说了!
25

RNA-seq流程需要进化啦!

Tophat 首次被发表已经是6年前

Cufflinks也是五年前的事情了

Star的比对速度是tophat的50倍,hisat更是star的1.2倍。

stringTie的组装速度是cufflinks的25倍,但是内存消耗却不到其一半。

Ballgown在差异分析方面比cuffdiff更高的特异性及准确性,且时间消耗不到cuffdiff的千分之一

Bowtie2+eXpress做质量控制优于tophat2+cufflinks和bowtie2+RSEM

Sailfish更是跳过了比对的步骤,直接进行kmer计数来做QC,特异性及准确性都还行,但是速度提高了25倍

kallisto同样不需要比对,速度比sailfish还要提高5倍!!!

参考:https://speakerdeck.com/stephenturner/rna-seq-qc-and-data-analysis-using-the-tuxedo-suite

24

用 GMAP/GSNAP软件进行RNA-seq的alignment

软件的解说ppt :http://www.mi.fu-berlin.de/wiki/pub/ABI/CompMethodsWS11/MHuska_GSNAP.pdf

软件的下载地址: http://research-pub.gene.com/gmap/
有研究者认为这个软件的比对效果要比tophat要好,虽然现在已经多出来了非常多的RNA-seq的alignment软件,我还是简单看看这个软件吧,它本来是2005就出来的一个专门比对低通量的est序列,叫GMAP,后来进化成了GSNAP
step1:下载安装GMAP/GSNAP
是一个标准的linux源码程序,安装之前一定要看readme  ,http://research-pub.gene.com/gmap/src/README
解压进去,然后源码安装三部曲,首先 ./configu  然后make 最后make install
会默认安装在 /usr/local/bin 下面,这里需要修改,因为你可能没有 /usr/local/bin 权限,安装到自己的目录,然后把它添加到环境变量!
step2 :准备数据
比对一般都只需要两个数据,一是索引好的参考基因组,另一个是需要比对的测序数据。
但是这个GSNAP,还需要对应的GTF注释文件。
首先需要参考基因组:虽然软件本身提供了一个hg19的参考基因组,并且已经索引好了Human genome, version hg19 (5.5 GB)(http://research-pub.gene.com/gmap/genomes/hg19.tar.gz) ,但是下载很慢,而且不是对所有版本的GSNAP都适用。所以我这里对我自己的参考基因组进行索引。
gmap_build -D ./ -d  my_hg19.fa
然后取ensemble下载hg19的gtf文件。
然后还需要把自己下载的gtf文件也构建索引,需要两个步骤
cat my_hg19.gtf |  ~/software/gmap-2011-10-16/util/gtf_splicesites > my_hg19.splicesites
cat  my_hg19.splicesites  |   iit_store -o my_hg19.gtf.index
然后拷贝需要比对的RNA-seq测序文件
step3: 运行程序
就是一步比对而已
gsnap
-D /home/jschnable/gsnap_indexes/
-d arabidopsisv10
--nthreads=50
-B 5
-s  /home/jschnable/gsnap_indexes/arabidopsisv10.iit
-n 2
-Q
--nofails
--format=sam temp.fastq
> results.sam
参数有点多,自己看看说明书吧http://qteller.com/RNAseq-analysis-recipe.pdf 讲的非常详细。
17

用DESeq进行差异分析的源代码

要保证当前文件夹下面有了742KO1.count等4个文件,就是用htseq对比对的bam文件进行处理后的输出文件

library(DESeq)
#加载数据
K1=read.table("742KO1.count",row.names=1)
K2=read.table("743KO2.count",row.names=1)
W1=read.table("740WT1.count",row.names=1)
W2=read.table("741WT2.count",row.names=1)
#列名
data=cbind(K1,K2,W1,W2)
#如果是htseq的结果,则删除data最后四行
n=nrow(data)
data=data

[c language="(-n+4:-n),"][/c]

#如果是bedtools的结果,取出统计个数列和行名
kk1=cbind(K1$V5)
rownames(kk1)=rownames(K1)
K1=kk1

#差异分析
colnames(data)=c("K1","K2","W1","W2")
type=rep(c("K","W"),c(2,2))
de=newCountDataSet(data,type)
de=estimateSizeFactors(de)
de=estimateDispersions(de)
res=nbinomTest(de,"K","W")

#res就是我们的表达量检验结果

到这里,理论上差异基因的分析已经结束啦!后面只是关于R的bioconductor包的一些简单结合使用而已

library(org.Mm.eg.db)

tmp=select(org.Mm.eg.db, keys=res$id, columns=c("ENTREZID","SYMBOL"), keytype="ENSEMBL")

#合并res和tmp
res=merge(tmp,res,by.x="ENSEMBL",by.y="id",all=TRUE)

#go
tmp=select(org.Mm.eg.db, keys=res$ENSEMBL, columns="GO", keytype="ENSEMBL")
ensembl_go=unlist(tapply(tmp[,2],as.factor(tmp[,1]),function(x) paste(x,collapse ="|"),simplify =F))

#为res加入go注释,
res$go=ensembl_go[res$ENSEMBL]#为res加入一列go

#写入all——data
all_res=res
write.csv(res,file="all_data.csv",row.names =F)

uniq=na.omit(res)#删除无效基因
sort_uniq=uniq[order(uniq$padj),]#按照矫正p值排序

#写入排序后的all_data
write.csv(res,file="all_data.csv",row.names =F)

#标记上下调基因
sort_uniq$up_down=ifelse(sort_uniq$baseMeanA>sort_uniq$baseMeanB,"up","down")
#交换上下调基因列位置
final_res=sort_uniq[,c(12,1:11)]
#写出最后数据
write.csv(final_res,file="final_annotation_gene_bedtools_sort_uniq.csv",row.names =F)

#然后挑选出padj值小于0.05的数据来做富集
tmp=select(org.Mm.eg.db, keys=sort_uniq[sort_uniq$padj<0.05,1], columns="ENTREZID", keytype="ENSEMBL")
diff_ENTREZID=tmp$ENTREZID
require(DOSE)
require(clusterProfiler)
diff_ENTREZID=na.omit(diff_ENTREZID)
ego <- enrichGO(gene=diff_ENTREZID,organism="mouse",ont="CC",pvalueCutoff=0.05,readable=TRUE)
ekk <- enrichKEGG(gene=diff_ENTREZID,organism="mouse",pvalueCutoff=0.05,readable=TRUE)
write.csv(summary(ekk),"KEGG-enrich.csv",row.names =F)
write.csv(summary(ego),"GO-enrich.csv",row.names =F)

 

21

RNA-seq流程对基因和转录本的表达量的计算

bedtools multicov和htseq-count都可以用来对基因和转录本的表达量的计算!!!

我们总共有四个样本,已经比对到小鼠的mm9基因组上面了,数据大小如下

RNA-seq流程对基因和转录本的表达量的计算111

然后对基因和转录本计数需要一些额外的信息,即各个基因及转录本的位置信息,gtf文件需要在UCSC等各大数据库下载

RNA-seq流程对基因和转录本的表达量的计算170

然后我们制作一个config文件配置我们的数据地址

cat sample_bam.config 可以看到文件内容如下

/data/mouse/ptan/740WT1.bam

/data/mouse/ptan/741WT2.bam

/data/mouse/ptan/742KO1.bam

/data/mouse/ptan/743KO2.bam

几个批处理文件名及内容分别如下

bedtools_multicov.sh  bedtools_multicov_transcript.sh  htseq.sh  htseq_transcript.sh

 

while read id

do

echo $id

new=`echo $id |cut -d"/" -f 5`

echo $new

bedtools multicov -bams $id -bed /data/mouse/mouse_mm9_gene.bed  > $new.gene.bedtools_multicov.count

done <$1

 

while read id

do

echo $id

new=`echo $id |cut -d"/" -f 5`

echo $new

bedtools multicov -bams $id -bed /data/mouse/mouse_mm9_transcript.bed  > $new.transcript.bedtools_multicov.count

done <$1

 

while read id

do

echo $id

new=`echo $id |cut -d"/" -f 5`

htseq-count -f bam $id /data/mouse/Mus_musculus.NCBIM37.67.gtf.chr  > $new.gene.htseq.count

done <$1

 

while read id

do

echo $id

new=`echo $id |cut -d"/" -f 5`

htseq-count -f bam --idattr transcript_id $id /data/mouse/Mus_musculus.NCBIM37.67.gtf.chr  > $new.transcript.htseq.count

done <$1

 

批量运行这些程序后就能对它们分别分情况进行计数,也能比较这两种计数方法的区别!

RNA-seq流程对基因和转录本的表达量的计算1201

 

可以看出区别还是很大的!!!

RNA-seq流程对基因和转录本的表达量的计算1219

我肯定没搞懂它们的原理,这完全就不一样,已经不是区别的问题了!!!

对于每个个体输出的计数文件,接下来就可以用DESeq等包来进行差异基因分析啦!

21

使用Bedtools对RNA-seq进行基因计数

以前是没有想过用这个软件的,直到有一个我的htseq无法对比对的bam文件进行基因计数(后来我才发现htseq无法计数的原因是gtf版本不同导致坐标不同,而且gtf对染色体编号没有加上chr),我简单搜索了一下,发现bedtools multicov也有类似的功能,所以我选择它来试试看!

首先注意它需要sort的bam文件及bam的index

bedtools multicov depends upon index BAM files in order to count the number of overlaps in each BAM file. As such, each BAM file should be position sorted (samtool sort aln.bam aln.sort) and indexed (samtools index aln.sort.bam) with either samtools or bamtools.

首先安装它:

wget https://github.com/arq5x/bedtools2/releases/download/v2.23.0/bedtools-2.23.0.tar.gz

解压开后

Make clean

Make all

然后就可以看到它的bin目录下全部是程序啦

Bedtools使用笔记639

命令很简单的

bedtools multicov [OPTIONS] -bams BAM1 BAM2 BAM3 ... BAMn -bed  <BED/GFF/VCF>

By default, multicov will report the count of alignments in each input BAM file that overlap.

例子:

bedtools multicov -bams aln1.bam aln2.bam aln3.bam -bed ivls-of-interest.bed

ivls-of-interest.bed这个文件是必须的,可能需要自己制作,其实用gtf文件也可以的

chr1 0   10000   ivl1

chr1 10000   20000   ivl2

chr1 20000   30000   ivl3

chr1 30000   40000   ivl4

输出结果前三列是坐标,第四列是基因名,跟我们的bed文件一样,只是最后三列是三个样本的计数,是添加上来的!

chr1 0       10000   ivl1    100 2234    0

chr1 10000   20000   ivl2    123 3245    1000

chr1 20000   30000   ivl3    213 2332    2034

chr1 30000   40000   ivl4    335 7654    0

 

同样是对gene的reads计数,bedtools的multicov工具与htseq-count的区别是

i'd guess it's due to the fact that htseq-count only reports one hit per aligned read assuming that read is aligned uniquely and does not overlap multiple features listed in your GTF. if an aligned read hits more than one feature in your GTF then it doesn't report that hit. bedtools gives you raw hits which includes every 1 hit for every intersection of every alignment with any features in the GTF no matter how many times it aligned or how many features it hit. you might think, "wow, htseq-count is dropping a lot of information". yes, it is! i've moved to using other tools to count hits to genes (RSEM/eXpress) since they disambiguate those ambiguous alignments and as a result you get counts for all of your aligned reads. in a genome with alternative splicing you lose too much data using htseq-count, in my opinion.

而且专门有个文献在讨论这个问题

http://barcwiki.wi.mit.edu/wiki/SOPs/rna-seq-diff-expressions

http://www.nature.com/nbt/journal/v31/n1/abs/nbt.2450.html

Differential analysis of gene regulation at transcript resolution with RNA-seq

下面我讲一个实际的例子

我的bam文件如下

Bedtools使用笔记2406

bedtools multicov -bams 740WT1.bam 741WT2.bam 742KO1.bam 743KO2.bam -bed mm9.bed

Bedtools使用笔记2491

得到的这个矩阵就可以去用DESeq包来进行差异分析啦!

18

R语言DESeq找差异基因

一:安装并加装该R包

安装就用source("http://bioconductor.org/biocLite.R") ;biocLite("DESeq")即可,如果安装失败,就需要自己下载源码包,然后安装R模块。

 

二.所需要数据

它的说明书指定了我们一个数据

source("http://bioconductor.org/biocLite.R") ;biocLite("pasilla")

安装了pasilla这个包之后,在这个包的安装目录就可以找到一个表格文件,就是我们的DESeq需要的文件。

C:\Program Files\R\R-3.2.0\library\pasilla\extdata\pasilla_gene_counts.tsv

说明书原话是这样的

The table cell in the i-th row and the j-th column of the table tells how many reads have been mapped to gene i in sample j.

一般我们需要用htseq-count这个程序对我们的每个样本的sam文件做处理计数,并合并这样的数据

下面这个是示例数据,第一列是基因ID号,后面的每一列都是一个样本。

图片1

de = newCountDataSet( pasillaCountTable, condition )  #根据我们的样本因子把基因计数表格读入成一个cds对象,这个newCountDataSet函数就是为了构建对象!

对我们构建好的de对象就可以直接开始找差异啦!非常简单的几步即可

de=estimateSizeFactors(de)

de=estimateDispersions(de)

res=nbinomTest(de,"K","W") #最重要的就是这个res表格啦!

uniq=na.omit(res)

我这里是对4个样本用htseq计数后的文件来做的,贴出完整代码吧

library(DESeq)

#首先读取htseq对bam或者sam比对文件的计数结果

K1=read.table("742KO1.count",row.names=1)

K2=read.table("743KO2.count",row.names=1)

W1=read.table("740WT1.count",row.names=1)

W2=read.table("741WT2.count",row.names=1)

data=cbind(K1,K2,W1,W2)

data=data[-c(43630:43634),]

#把我们的多个样本计数结果合并起来成数据框,列是不同样本,行是不同基因

colnames(data)=c("K1","K2","W1","W2")

type=rep(c("K","W"),c(2,2))

#构造成DESeq的对象,并对分组样本进行基因表达量检验

de=newCountDataSet(data,type)

de=estimateSizeFactors(de)

de=estimateDispersions(de)

res=nbinomTest(de,"K","W")

#res就是我们的表达量检验结果

library(org.Mm.eg.db)

tmp=select(org.Mm.eg.db, keys=res$id, columns="GO", keytype="ENSEMBL")

ensembl_go=unlist(tapply(tmp[,2],as.factor(tmp[,1]),function(x) paste(x,collapse ="|"),simplify =F))

#首先输出所有的计数数据,加上go注释信息

all_res=res

res$go=ensembl_go[res$id]

write.csv(res,file="all_data.csv",row.names =F)

#然后输出有意义的数据,即剔除那些没有检测到表达的基因

uniq=na.omit(res)

sort_uniq=uniq[order(uniq$padj),]

write.csv(sort_uniq,file="sort_uniq.csv",row.names =F)

#然后挑选出padj值小于0.05的差异基因数据来做富集,富集用的YGC的两个包,在我前面的博客已经详细说明了!

tmp=select(org.Mm.eg.db, keys=sort_uniq[sort_uniq$padj<0.05,1], columns="ENTREZID", keytype="ENSEMBL")

diff_ENTREZID=tmp$ENTREZID

require(DOSE)

require(clusterProfiler)

diff_ENTREZID=na.omit(diff_ENTREZID)

ego <- enrichGO(gene=diff_ENTREZID,organism="mouse",ont="CC",pvalueCutoff=0.01,readable=TRUE)

ekk <- enrichKEGG(gene=diff_ENTREZID,organism="mouse",pvalueCutoff=0.01,readable=TRUE)

write.csv(summary(ekk),"KEGG-enrich.csv",row.names =F)

write.csv(summary(ego),"GO-enrich.csv",row.names =F)

 

10

RNA-seq比对软件HISAT说明书

取代bowtie+tophat进行RNA-seq比对

HISAT全称为Hierarchical Indexing for Spliced Alignment of Transcripts,由约翰霍普金斯大学开发。它取代Bowtie/TopHat程序,能够将RNA-Seq的读取与基因组进行快速比对。这项成果发表在3月9日的《Nature Methods》上。

HISAT利用大量FM索引,以覆盖整个基因组。以人类基因组为例,它需要48,000个索引,每个索引代表~64,000 bp的基因组区域。这些小的索引结合几种比对策略,实现了RNA-Seq读取的高效比对,特别是那些跨越多个外显子的读取。尽管它利用大量索引,但HISAT只需要4.3 GB的内存。这种应用程序支持任何规模的基因组,包括那些超过40亿个碱基的。

HISAT软件可从以下地址获取:http://ccb.jhu.edu/software/hisat/index.shtml。

首先,我们安装这个软件!

Wget http://ccb.jhu.edu/software/hisat/downloads/hisat-0.1.5-beta-source.zip

官网下载的是源码包,需要make一下,make之后目录下面就多了很多程序,绿色的那些都是,看起来是不是很眼熟呀!!!

哈哈,这完全就是bowtie的模拟版本!!!

HISAT取代bowtie+tophat进行RNA-seq比对1222

也可以从github里面下载,wget https://codeload.github.com/infphilo/hisat/zip/master

下载后直接解压即可使用啦。当然这个软件本身也有着详尽的说明书

http://ccb.jhu.edu/software/hisat/manual.shtml

然后就是准备数据,它跟tophat一样的功能。就是把用RNA-seq方法测序得到的fastq文件比对到参考基因组上面,所以就准这两个文件了哦

接下来是运行程序!

说明书上面写着分成两个步骤,构建索引和比对。

这个软件包模仿bowtie自带了一个example数据,而且它的说明书也是针对于那个example来的,我也简单运行一下。

$HISAT_HOME/hisat-build $HISAT_HOME/example/reference/22_20-21M.fa 22_20-21M_hisat

构建索引的命令如上,跟bowtie一样我修改了一下

/home/jmzeng/hoston/RNA-soft/hisat-0.1.5-beta/hisat-build 22_20-21M.fa  my_hisat_index

连日志都跟bowtie一模一样,哈哈,可以看到我们的这个参考fasta文件 22_20-21M.fa 就变成索引文件啦,索引还是很多的!

HISAT取代bowtie+tophat进行RNA-seq比对1871

然后就是比对咯,还是跟bowtie一样

$HISAT_HOME/hisat -x 22_20-21M_hisat -U $HISAT_HOME/example/reads/reads_1.fq -S eg1.sam

我的命令是

/home/jmzeng/hoston/RNA-soft/hisat-0.1.5-beta/hisat -x  my_hisat_index -U ../reads/reads_1.fq  -S reads1.sam

1000 reads; of these:

1000 (100.00%) were unpaired; of these:

0 (0.00%) aligned 0 times

1000 (100.00%) aligned exactly 1 time

0 (0.00%) aligned >1 times

100.00% overall alignment rate

哈哈,到这里。这个软件就运行完毕啦!!!是不是非常简单,只有你会用bowtie,这个就没有问题。当然啦,软件还是有很多细节是需要调整的。我下面就简单讲一个实际的例子哈!

首先,我用了1.5小时把4.6G的小鼠基因组构建了索引

/home/jmzeng/hoston/RNA-soft/hisat-0.1.5-beta/hisat-build  Mus_musculus.GRCm38.fa.fa mouse_hisat_index

HISAT取代bowtie+tophat进行RNA-seq比对2512

然后对我的四个测序文件进行比对。

for i in *fq

do

/home/jmzeng/hoston/RNA-soft/hisat-0.1.5-beta/hisat  -x  /home/jmzeng/hoston/mouse/mouse_hisat_index  \

-p 30 -U  $i.trimmed.single  -S ./hisat_out/${i%.*}.sam

done

它运行的速度的确要比tophat快好多,太可怕的速度!!!!至于是否多消耗了内存我就没有看了

4.6G的小鼠,5G的测序数据,我只用了五个核,居然十分钟就跑完了!

然后听群友说是因为没有加 --known-splicesite-infile <path>这个参数的原因,没有用gtf文件来指导我们的RNA数据的比对,这样是不对的!

需要用下面这个脚本把gtf文件处理一下,然后导入什么那个参数来指导RNA比对。

extract_splice_sites.py genes.gtf > splicesites.txt

但是我报错了,错误很奇怪,没解决,但是我换了个 extract_splice_sites.py  程序,就可以运行啦!之前是HISAT 0.1.5-beta release 2/25/2015里面的python程序,后来我换做了github里面的就可以啦!

/home/jmzeng/hoston/RNA-soft/hisat-master/extract_splice_sites.py Mus_musculus.GRCm38.79.gtf >mouse_splicesites.txt

HISAT取代bowtie+tophat进行RNA-seq比对3218

21192819 reads; of these:
21192819 (100.00%) were unpaired; of these:
14236834 (67.18%) aligned 0 times
5437800 (25.66%) aligned exactly 1 time
1518185 (7.16%) aligned >1 times

感觉没有变化,不知道为什么?

21192819 reads; of these:

21192819 (100.00%) were unpaired; of these:

14236838 (67.18%) aligned 0 times

5437793 (25.66%) aligned exactly 1 time

1518188 (7.16%) aligned >1 times

32.82% overall alignment rate

发表这个软件的文献本身也把这个软件跟其它软件做了详尽的对比

http://www.nature.com/nmeth/journal/v12/n4/full/nmeth.3317.html

Program Run time (min) Memory usage (GB)
Run times and memory usage for HISAT and other spliced aligners to align 109 million 101-bp RNA-seq reads from a lung fibroblast data set. We used three CPU cores to run the programs on a Mac Pro with a 3.7 GHz Quad-Core Intel Xeon E5 processor and 64 GB of RAM.
HISATx1 22.7 4.3
HISATx2 47.7 4.3
HISAT 26.7 4.3
STAR 25 28
STARx2 50.5 28
GSNAP 291.9 20.2
OLego 989.5 3.7
TopHat2 1,170 4.3

 

 

 

参考:http://www.plob.org/2015/03/20/8980.html

http://nextgenseek.com/2015/03/hisat-a-fast-and-memory-lean-rna-seq-aligner/

 

10

RNA-seq的比对软件star说明书

类似于tophat的软件

首先当然是下载软件啦!

两个地方可以下载,一个是谷歌code中心,被墙啦,另一个是github,我的最爱。

wget https://codeload.github.com/alexdobin/STAR/zip/master

2013-RNA-seq的star-类似于tophat的软件说明书217

解压即可使用啦,其中程序在bin目录下面,根据自己的平台调用即可!

然后doc里面还有个pdf的说明文档,写的非常清楚,我也是看着那个文档学的这个软件!

接下来就是准备数据啦!

既然是类似于tophat一样的比对软件,当然是准备参考基因组和测序数据咯,毫无悬念。

然后 该软件也给出了一些测试数据

ftp://ftp2.cshl.edu/gingeraslab/tracks/STARrelease/2.1.4/

然后就是运行程序的命令!

分为两个步骤:首先构建索引,然后比对即可,中间的参数根据具体需要可以细调!

构建索引时候,软件说明书给的例子是

The basic options to generate genome indices are as follows:
--runThreadN NumberOfThreads
--runMode genomeGenerate
--genomeDir /path/to/genomeDir
--genomeFastaFiles /path/to/genome/fasta1 /path/to/genome/fasta2 ...
--sjdbGTFfile /path/to/annotations.gtf
--sjdbOverhang ReadLength-1

我模仿了一下。对我从ensembl ftp里面下载的老鼠基因组构建了索引

/home/jmzeng/hoston/RNA-soft/STAR-master/bin/Linux_x86_64/STAR  \

--runThreadN 30 #我的服务器还比较大,可以使用30个CPU  \

--runMode genomeGenerate \

--genomeDir  /home/jmzeng/hoston/mouse/STAR-mouse #构建好的索引放在这个目录 \

--genomeFastaFiles  /home/jmzeng/hoston/mouse/Mus_musculus.GRCm38.fa.fa \

--sjdbGTFfile /home/jmzeng/hoston/mouse/Mus_musculus.GRCm38.79.gtf \

--sjdbOverhang 284   #我的测序数据长短不一,最长的是285bp

当然注释的地方你要删除掉才行呀,因为cpu用的比较多。

2013-RNA-seq的star-类似于tophat的软件说明书1302

算一算时间,对4.6G的小鼠基因组来说,半个小时算是非常快的了!Bowtie2的index要搞两个多小时。

然后就是比对咯。这也是很简单的,软件说明书给的例子是

The basic options to run a mapping job are as follows:
--runThreadN NumberOfThreads
--genomeDir /path/to/genomeDir
--readFilesIn /path/to/read1 [/path/to/read2]

我稍微理解了一下参数,然后写出了自己的命令。

fq=740WT1.fq.trimmed.single

mkdir  740WT1_star

/home/jmzeng/hoston/RNA-soft/STAR-master/bin/Linux_x86_64/STAR  \

--runThreadN 20  \

--genomeDir  /home/jmzeng/hoston/mouse/STAR-mouse   \

--readFilesIn  $fq \

--outFileNamePrefix  ./740WT1_star/740WT1

如果输出文件需要被cufflinks套装软件继续使用。就需要用一下参数

Cufflinks/Cuffdiff require spliced alignments with XS strand attribute, which STAR will generate with --outSAMstrandField intronMotif option.

还有--outSAMtype参数可以修改输出比对文件格式,可以是sam也可以是bam,可以是sort好的,也可以是不sort的。

最后是输出文件解读咯!

其实没什么好解读的,输出反正就是sam类似的比对文件咯,如果还有其它文件,需要自己好好解读说明书啦。我就不废话了!

 

值得一提的是,该程序提供了2次map的建议

The basic idea is to run 1st pass of STAR mapping with the usual parameters , then collect the junctions detected in the first pass, and use them as ”annotated” junctions for the 2nd pass mapping.

在对RNA-seq做snp-calling的时候可以用到,尤其是GATK官方还给出了教程,大家可以好好学习学习!

http://www.broadinstitute.org/gatk/guide/article?id=3891