Tag Archives: RNA

九 25

RNA-seq流程需要进化啦！

Posted on by

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

Posted on by

软件的解说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 讲的非常详细。
五 05

RNA-seq完整学习手册！

Posted on by

需耗时两个月！里面网盘资料如果过期了，请直接联系我1227278128，或者我的群201161227，所有的资源都可以在 http://pan.baidu.com/s/1jIvwRD8 此处找到

搜索可以得到非常多的流程，我这里简单分享一些，我以前搜索到的文献。

 

RNA-seq完整学习手册141

北大也有讲RNA-seq的原理

链接：http://pan.baidu.com/s/1kTmWmv9 密码：6yaz

甚至，我还有个华大的培训课程！！！这可是5天的培训教程哦，好像当初还花了五千多块钱的资料！！！

链接：http://pan.baidu.com/s/1nt5OV5B 密码：gyul

RNA-seq完整学习手册294

优酷也有视频，可以自己搜索看看

RNA-seq完整学习手册312

然后还有几个pipeline，就是生信的分析流程，即使你啥都不会，按照pipeline来也不是问题啦

export PATH=/share/software/bin:$PATH

bowtie2-build ./data/GRCh37_chr21.fa  chr21

tophat -p 1 -G ./data/genes.gtf -o P460.thout chr21 ./data/P460_R1.fq  ./data/P460_R2.fq

tophat -p 1 -G ./data/genes.gtf -o C460.thout chr21 ./data/C460_R1.fq  ./data/C460_R2.fq

cufflinks -p 1 -o P460.clout P460.thout/accepted_hits.bam

cufflinks -p 1 -o C460.clout C460.thout/accepted_hits.bam

samtools  view  -h  P460.thout/accepted_hits.bam  >  P460.thout/accepted_hits.sam

samtools  view  -h  C460.thout/accepted_hits.bam  >  C460.thout/accepted_hits.sam

echo ./P460.clout/transcripts.gtf > assemblies.txt

echo ./C460.clout/transcripts.gtf >> assemblies.txt

cuffmerge -p 1 -g ./data/genes.gtf -s ./data/GRCh37_chr21.fa  assemblies.txt

cuffdiff -p 1 -u merged_asm/merged.gtf  -b ./data/GRCh37_chr21.fa  -L P460,C460 -o P460-C460.diffout P460.thout/accepted_hits.bam C460.thout/accepted_hits.bam

samtools  index  P460.thout/accepted_hits.bam

samtools  index  C460.thout/accepted_hits.bam

 

和另外一个

#!/bin/bash

# Approx 75-80m to complete as a script

cd ~/RNA-seq

ls -l data

 

tophat --help

 

head -n 20 data/2cells_1.fastq

 

time tophat --solexa-quals \

-g 2 \

--library-type fr-unstranded \

-j annotation/Danio_rerio.Zv9.66.spliceSites\

-o tophat/ZV9_2cells \

genome/ZV9 \

data/2cells_1.fastq data/2cells_2.fastq                  # 17m30s

 

time tophat --solexa-quals \

-g 2 \

--library-type fr-unstranded \

-j annotation/Danio_rerio.Zv9.66.spliceSites\

-o tophat/ZV9_6h \

genome/ZV9 \

data/6h_1.fastq data/6h_2.fastq                          # 17m30s

 

samtools index tophat/ZV9_2cells/accepted_hits.bam

samtools index tophat/ZV9_6h/accepted_hits.bam

 

cufflinks --help

time cufflinks  -o cufflinks/ZV9_2cells_gff \

-G annotation/Danio_rerio.Zv9.66.gtf \

-b genome/Danio_rerio.Zv9.66.dna.fa \

-u \

--library-type fr-unstranded \

tophat/ZV9_2cells/accepted_hits.bam                  # 2m

 

 

time cufflinks  -o cufflinks/ZV9_6h_gff \

-G annotation/Danio_rerio.Zv9.66.gtf \

-b genome/Danio_rerio.Zv9.66.dna.fa \

-u \

--library-type fr-unstranded \

tophat/ZV9_6h/accepted_hits.bam                      # 2m

 

# guided assembly

time cufflinks  -o cufflinks/ZV9_2cells \

-g annotation/Danio_rerio.Zv9.66.gtf \

-b genome/Danio_rerio.Zv9.66.dna.fa \

-u \

--library-type fr-unstranded \

tophat/ZV9_2cells/accepted_hits.bam                  # 16m

 

 

time cufflinks  -o cufflinks/ZV9_6h \

-g annotation/Danio_rerio.Zv9.66.gtf \

-b genome/Danio_rerio.Zv9.66.dna.fa \

-u \

--library-type fr-unstranded \

tophat/ZV9_6h/accepted_hits.bam                      # 13m

 

 

time cuffdiff -o cuffdiff/ \

-L ZV9_2cells,ZV9_6h \

-T \

-b genome/Danio_rerio.Zv9.66.dna.fa \

-u \

--library-type fr-unstranded \

annotation/Danio_rerio.Zv9.66.gtf \

tophat/ZV9_2cells/accepted_hits.bam \

tophat/ZV9_6h/accepted_hits.bam                        # 7m

 

head -n 20 cuffdiff/gene_exp.diff

 

sort -t$'\t' -g -k 13 cuffdiff/gene_exp.diff \

> cuffdiff/gene_exp_qval.sorted.diff

 

head -n 20 cuffdiff/gene_exp_qval.sorted.diff

三 19

转录组-TransDecoder-对trinity结果进行注释

Posted on by

   一：下载安装该软件

下载安装该软件：  wget https://codeload.github.com/TransDecoder/TransDecoder/tar.gz/2.0.1

解压进入该目录，查看里面的文件

make一下就可以用了，看起来好像是依赖于perl模块的

转录组-TransDecoder-预测ORF420

这个TransDecoder.LongOrfs就是我们这次需要的程序，查看该程序，的确真是一个perl程序，看来perl还是蛮有用的。

二：准备数据

它里面有个测试数据，是比较全面的，也比较复杂，我就不贴出来了，反正我是那trinity组装好的fasta格式的转录组数据来预测ORF的。

三：运行命令

它给的测试命令也很复杂

## generate alignment gff3 formatted output

../util/cufflinks_gtf_to_alignment_gff3.pl transcripts.gtf > transcripts.gff3

 

## generate transcripts fasta file

../util/cufflinks_gtf_genome_to_cdna_fasta.pl transcripts.gtf test.genome.fasta > transcripts.fasta 

 

## Extract the long ORFs

../TransDecoder.LongOrfs -t transcripts.fasta

当然我们只需要看最后一步，这是重点

我这里是直接对我们的trinity组装好的转录本进行预测ORF

/home/jmzeng/bio-soft/TransDecoder/TransDecoder.LongOrfs  -t Trinity.fasta

命令很简单

转录组-TransDecoder-预测ORF1471

输出来的文件就有预测的蛋白文件，这个文件是trinotate对转录本进行注释所必须的文件

转录组-TransDecoder-预测ORF1714

 

四：输出文件解读

longest_orfs.cds  这个是预测到的cds碱基序列，

longest_orfs.gff3  这个是预测得到的gff文件

longest_orfs.pep   这个就是预测得到的蛋白文件