高通量RNA测序已经能够更高效地检测融合转录本，但是融合检测的技术和相关软件通常产生高错误发现率。而一个自动整合RNA数据和已知基因组特征的可视化框架对于结果的检验是有帮助的。2017年发布的一个bioconductor包，chimeraviz就可以做到自动创建嵌合RNA可视化。
支持来自9种不同融合发现工具（deFuse、EricScript、InFusion、JAFFA、FusionCatcher、FusionMap、PRADA、SOAPfuse和STAR-FUSION）的输入。 Continue reading
该软件早在2016年就公布了，发表在biorxiv预印本上面，但直到2017年的双11，才发表在NG上面，文章是 : Annotation-free quantification of RNA splicing using LeafCutter 最大的特点应该是不需要参考基因组的基因注释信息了吧，就是gtf/gff文件可以省略，当然，比对还是需要的。它还有另外一个非常重要的功能，splicing quantitative trait loci (sQTLs) 但是跟我目前关系不大， 就不介绍了。 Continue reading
可变剪切是指mRNA前体以多种方式将exon连接在一起的过程。 由于可变剪切使一个基因产生多个mRNA转录本，不同mRNA可能翻译成不同蛋白。
转录组一般是指从细胞或组织的基因组所转录出来的RNA的总和，包括编码蛋白质的mRNA和各种非编码RNA（rRNA,tRNA,snRNA,snoRNA,lncRNA,microRNA等）。真核生物的基因结构是不连续的，如下图：
其实这个植物是拟南芥，所以跟人类研究的数据处理大同小异。
转录组测序的研究对象为特定细胞在某一功能状态下所能转录出来的所有 RNA 的总和，包括 mRNA 和非编码 RNA 。通过转录组测序，能够全面获得物种特定组织或器官的转录本信息，从而进行转录本结构研究、变异研究、基因表达水平研究以及全新转录本发现等研究。 Continue reading
早在去年九月，我就写个博文说 RNA-seq流程需要进化啦！ http://www.bio-info-trainee.com/1022.html ，主要就是进化成hisat2+stringtie+ballgown的流程，但是我一直没有系统性的讲这个流程，因为我觉真心木有用。我只用了里面的hisat来做比对而已！但是群里的小伙伴问得特别多，我还是勉为其难的写一个教程吧，你们之间拷贝我的代码就可以安装这些软件的！然后自己找一个测试数据，我的脚本很容易用的！ Continue reading
我以前写过bedtools和htseq-counts的教程，它们都可以用来对比对好的bam文件进行计数，正好群里有小伙伴问我它们的区别，我就简单做了一个比较，大家可以先看看我以前写的软件教程。写的有的挫：
言归正传，我这里没精力去探究它们的具体原理，只是看看它们数一个read是否属于某个基因的时候，区别在哪里，大家看下图： Continue reading
我这里拿的是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细胞系，处理详情如下：
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
下载的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
接下来所有的分析就基于此数据啦
有网友咨询过对于没有参考基因组或者转录组的物种，如何做RNA-seq分析。我觉得这个问题太大了，而且我还真的对这个没有经验。但是我以前看到过一篇文献，里面提到过一个非常全面的转录组 de novo组装注释流程，所以我摘抄了文章里面的生物信息学处理部分，分享给大家： Continue reading
开发单位：华大，SOAP系列软件套装！
二，输入数据准备
6.5K Jun 15 2009 cytoBand.txt.gz3.0G Oct 12 2012 hg19.fa2.5M Mar 15 10:30 HGNC_Gene_Family_dataset38M Feb 8 2014 Homo_sapiens.GRCh37.75.gtf.gz202 Jan 19 16:07 HumanRef_refseg_symbols_relationship.list
文件下载地址，作者已经给出了！
这一步耗时很长，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.fagtf_annotation_file: /home/jmzeng/biosoft/SOAPfuse/db_v1.27/raw/Homo_sapiens.GRCh37.75.gtf.gzChr_Bandregion_file: /home/jmzeng/biosoft/SOAPfuse/db_v1.27/raw/cytoBand.txt.gzHGNC_gene_family_file: /home/jmzeng/biosoft/SOAPfuse/db_v1.27/raw/HGNC_Gene_Family_datasetgtf_segname2refseg_list: /home/jmzeng/biosoft/SOAPfuse/db_v1.27/raw/HumanRef_refseg_symbols_relationship.list
DB_db_dir = /DATABASE_DIR/PG_pg_dir = /TOOL_DIR/source/binPS_ps_dir = /TOOL_DIR/sourcePD_all_out = /out_directory/PA_all_fq_postfix = PostFix
三，运行命令
perl SOAPfuse-RUN.pl -c <config_file> -fd <WHOLE_SEQ-DATA_DIR> -l <sample_list> -o <out_directory> [Options]
运行的非常慢！！！
四，数据结果解读
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
软件的解说ppt ：http://www.mi.fu-berlin.de/wiki/pub/ABI/CompMethodsWS11/MHuska_GSNAP.pdf
要保证当前文件夹下面有了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)
bedtools multicov和htseq-count都可以用来对基因和转录本的表达量的计算！！！
我们总共有四个样本，已经比对到小鼠的mm9基因组上面了，数据大小如下
然后对基因和转录本计数需要一些额外的信息，即各个基因及转录本的位置信息，gtf文件需要在UCSC等各大数据库下载
然后我们制作一个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
批量运行这些程序后就能对它们分别分情况进行计数，也能比较这两种计数方法的区别！
可以看出区别还是很大的！！！
我肯定没搞懂它们的原理，这完全就不一样，已经不是区别的问题了！！！
对于每个个体输出的计数文件，接下来就可以用DESeq等包来进行差异基因分析啦！
以前是没有想过用这个软件的，直到有一个我的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 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 multicov -bams 740WT1.bam 741WT2.bam 742KO1.bam 743KO2.bam -bed mm9.bed
得到的这个矩阵就可以去用DESeq包来进行差异分析啦！
安装就用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号，后面的每一列都是一个样本。
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)
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的模拟版本！！！
也可以从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 就变成索引文件啦，索引还是很多的！
然后就是比对咯，还是跟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
然后对我的四个测序文件进行比对。
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
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/
首先当然是下载软件啦！
两个地方可以下载，一个是谷歌code中心，被墙啦，另一个是github，我的最爱。
wget https://codeload.github.com/alexdobin/STAR/zip/master
解压即可使用啦，其中程序在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用的比较多。
算一算时间，对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