十二 09

【直播】我的基因组(十一):测序数据的比对

上一次直播中,我们对拿到手的测序数据进行了质控,测序数据的质量已经得到了保证。那么接下来就可以把它拿来与参考基因组比对了,这里我们先用参考基因组hg19,大家可以参照【直播】我的基因组(五):测试数据及参考基因组的准备来下载参考基因组hg19,我这里选择的是UCSC提供的hg19。然后安装bwa软件进行比对,可以参考【直播】我的基因组(四):计算资源的准备来安装,以及对hg19建立索引。

我们首先简单讲一下为什么要进行比对以及比对过程是怎样的?

可以看到我们到手的测序数据格式是fastq,每条reads都是150个碱基组成,如果只看这fastq,我们没办法得知它的意义,参考基因组那么大(人类约30亿个碱基),这个reads在我们基因组的哪里呢?

简单解释一下,假设人类基因组是123456789,如果我们的测序得到的reads是123,那么我们很明显知道这条reads在基因组的第一个位置,跨越了3个长度,如果我们的reads是567,那么我们也可以看到它在基因组的第5个位置。如果我们看到了一个reads是235567,同样我们也很容易看到它在基因组第2位置,但是在第3个位置参考是4,它却是5,这里可能是测序错误,也可能是这个reads真的变异了!

但是我们的参考基因组远远没有那么简单,30亿个碱基的庞大数目,测序的一条reads也有150个碱基,仅仅用肉眼观察基本不可能判断出它到底在哪里。但并非一定观察不到,如果你有多的不可计的时间及精力的话,手工比对穷极一生来搞定一条reads的比对就很不容易了(当然肯定不会有人这么傻,这里只是说数据量真的很大而已)。然而在我们手上可是有8.9亿条reads,所以我们需要借助计算机来进行比对,现在比较流行的基因组比对工具是bwa和bowtie,它俩的算法不一样,但是我们不需要了解那么具体,只需要知道它可以把我们的fastq测序文件通过与参考基因组的比对生成sam格式(自行搜索了解该格式)的比对结果文件(如下),从sam文件中,我们可以看到每条reads在参考基因组的位置,这条reads是在哪一条染色体,又是在这条染色体的哪个位置就可以一目了然。

对于比对的结果,我们可以用IGV可视化查看,还可以手动查看每个基因的比对情况:

下面我简单讲一下代码

1,下载hg19基因组

cd ~/reference

mkdir -p genome/hg19  && cd genome/hg19

nohup wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz &

tar zvfx chromFa.tar.gz

cat *.fa > hg19.fa

rm chr*.fa

首先要理解linux基础命令,在我们的服务器上面新建好目录,找到hg19的下载链接,用linux自带的wget下载,因为文件太大,所以我们用nohup放在后台下载。下载后是压缩文件 chromFa.tar.gz,在linux里面需要用tar zvfx 来解压tar.gz文件即可。解压开后是一个个文件,需要用cat合并!最终效果如下:

2,安装bwa软件

## Download and install BWA

cd ~/biosoft

mkdir bwa &&  cd bwa

#http://sourceforge.net/projects/bio-bwa/files/

wget https://sourceforge.net/projects/bio-bwa/files/bwa-0.7.15.tar.bz2

tar xvfj bwa-0.7.15.tar.bz2 # x extracts, v is verbose (details of what it is doing), f skips prompting for each individual file, and j tells it to unzip .bz2 files

cd bwa-0.7.15

make

~/biosoft/bwa/bwa-0.7.15/bwa

我所有的软件都安装在自己的home目录下面的biosoft文件夹。同样,也是找的bwa的下载地址,然后解压,然后直接make即可。很多人的服务器会报错zlib.h缺少的问题,看我以前的教程:http://www.bio-info-trainee.com/518.html ,缺少什么你就安装什么,但是缺少的东西需要安装到系统环境变量,但是我的bwa是直接安装到自己的目录,所以我用全路径在调用该软件。如果你的这个命令~/biosoft/bwa/bwa-0.7.15/bwa  能够显示下面的help文档,说明你已经安装成功啦~

3,对hg19参考基因组用bwa构建索引

cd ~/reference

mkdir -p index/bwa && cd index/bwa

nohup time ~/biosoft/bwa/bwa-0.7.15/bwa index   -a bwtsw   -p ~/reference/index/bwa/hg19  ~/reference/genome/hg19/hg19.fa 1>hg19.bwa_index.log 2>&1   &

代码很简单,就是新建好一个文件夹来存放我们的参考基因组的索引,我这里选择的是我的home目录下面的reference/index/bwa/ 文件夹,可以看到如下内容:

我还是用了nohup把这个命令挂在后台,防止掉线,因为要运行2个小时左右,我加上time命令可以看到运行时间,我用了bwa的index模式来索引参考基因在,具体bwa用法可以自己看文档,但是我们只需要学会索引及比对就好了。有点类似于window下面的软件有一个个菜单栏一样,需要自己的鼠标点击来实现一个个功能,在linux下面就是把命令准备好,然后运行。

4.把fastq文件比对到参考基因组

for i in $(seq 1 6) ;do (nohup ~/biosoft/bwa/bwa-0.7.15/bwa  mem -t 5 -M ~/reference/index/bwa/hg19  KPGP-00001_L${i}_R1.fq.gz KPGP-00001_L${i}_R2.fq.gz 1>KPGP-00001_L${i}.sam 2>KPGP-00001_L${i}.bwa.align.log &);done

这个命令就一句话,但是里面的信息量非常大, 需要熟练掌握linux命令以及shell脚本的语法,但是解析起来也很简单,就是因为我们的fastq文件命名是有规律的,根据规律我构造出一个循环命令,里面的i这个变量会自动扩展成1,2,3,4,5,6依次来用bwa  mem 模式来比对,因为是PE150测序,所以选择这个模式,-M就是选择我们上一步构建好的参考基因组,最后面的 1> 和2>是把软件运行结果输出来,分别是标准输出和标准错误输出,大家可以自行搜索。如果fastq文件的命名发生变化,这个shell脚本是运行不了的,需要临时构建,自己得掌握脚本编写,不然就一个个的比对,手动。

大家可以去看【直播】我的基因组(七):从整体理解全基因组测序数据的变异位点,来了解这个命令的运行结果。

请扫描以下二维码关注我们,获取直播系列的所有帖子!

1

十一 23

关于multiple mapping我想说的

很多时候,我们都要选取unique mapped的reads,尤其是在RNA-seq和CHIP-seq的时候,但是如何保留,各种教程都不一致,我稍微总结了一下,是因为使用的比对工具不一样导致的!但是主要都反应在sam文件的一系列tag里面~

首先对bwa来说,如果它遇到一个reads可以比对到参考基因在的多个序列,只会随机的选取一个位置来输出到sam文件,但是会加上一个tag是XS:I:<N>来告诉我们第二好的比对情况的比对得分是多少,bowtie也是一样。但是它们都有参数来决定是否只对每个reads输出一条信息,还是输出全部的信息,在bwa是-a的参数,在bowtie里面是-m参数。

但是bowtie2里面取消了这个参数,它们都必须用XS:I:<N>这个tag来挑选unique mapped的reads

但是如果是用hisat来比对的话,决定是否是唯一比对的却是NH这个tag信息。默认情况下一条reads可以输出多条比对结果。

我想起了再补充吧,其实应该找几个例子用IGV看看,就明白了,可是我暂时没有时间了,只是觉得这个很重要,就提一下。

 

04

用 SHRiMP 来比对color space的数据

无意中接触了AB 5500xl Genetic Analyzer这个测序仪的数据,就是传说中的solid格式,也就是color space的测序数据 ,虽然拿到的测序数据也是fastq格式的, 4行代表一条read,但是第二行已经不是在是碱基序列啦,而是color的编码。Colors may be encoded either as numbers (0=blue, 1=green, 2=orange, 3=red) or as characters A/C/G/T (A=blue, C=green, G=orange, T=red).我们通常称为csfastq格式。
对于这种数据的处理,一般的比对软件是hold 不住的,我查了一下,SHRiMP,sequel和BFAST ,bowtie,是可以处理这种csfastq格式数据的比对的, 我这里简单使用了最出名的SHRiMP 。

Continue reading

05

华大soap系列的比对软件

也不知道是什么原因,对国产软件总是提不起兴趣,所以尽管SOAP系列都已经发展到了十几个软件了,我依然没有去试用一下。

软件下载:
官网直接找到:http://soap.genomics.org.cn/
SOAPaligner/soap2 is a member of the SOAP (Short Oligonucleotide Analysis Package).
很久以前,大家说soap其实指的是类似于bwa这样的比对工具,但是后来这个工具箱丰富了,所以我们现在如果只看比对工具,要看的是SOAPaligner
我是linux系统,用wget下载:wget http://soap.genomics.org.cn/down/soap2.21release.tar.gz
解压,由于下载是可执行程序,就不需要安装啦!
1
安装之后把该软件添加到环境变量!
输入数据:
这里选择两个网络上的测试数据:
如果是真想用这个软件的话,需要参考基因组和测序数据,这个链接貌似已经年久失修啦~!
# download a test reference genome (TAIR9 Chromosome 1)
wget 
http://biocluster.ucr.edu/~tbackman/query.fastq 
# download some test Illumina reads from Arabidopsis

运行命令:

2bwt-builder genome.fasta
   # create binary of reference genome
soap -a query.fastq -D genome.fasta.index -o output.soap
   # align query to genome and store output

结果解读:

由于测试数据没有下载下来,我安装了软件就懒得玩了,其实正经的来讲,应该写一个详细的测评,包括软件运行速度,比对准确率,等等,不过那样做就是发paper的节奏了,我随便玩玩,就算啦。
不过soap是一直在更新的,所以我相信他比对的结果,肯定是sam格式的。
所以结果就不用解读啦!
05

很老的比对软件maq

MAQ在2008年还是蛮火的,但是现在基本都是BWA和bowtie的天下了。
就当怀念一下它吧,给它写一个教程!
软件下载:
官网直接找到:http://maq.sourceforge.net/
解压,很容易观察到是C++源码,所以用源码安装三部曲来安装
tar jxvf software.tar.bz2
cd software
./configure --prefix=$path
make
make test
安装之后把该软件添加到环境变量!
输入数据:
这里选择两个网络上的测试数据:
如果是真想用这个软件的话,需要参考基因组和测序数据,这个链接貌似已经年久失修啦~!
# download a test reference genome (TAIR9 Chromosome 1)
wget 
http://biocluster.ucr.edu/~tbackman/query.fastq 
# download some test Illumina reads from Arabidopsis

运行命令:

maq # inspect command line options
maq fasta2bfa genome.fasta genome.bfa
   # create binary of reference genome
maq fastq2bfq query.fastq readBinary.bfq
   # create a binary of dataset
maq match out.map genome.bfa readBinary.bfq
# align query to genome and store output

结果解读:

我在想,这个MAQ软件发明之前,好像还没有SAM文件格式的定义,那么它的结果out.map肯定不是sam格式的。
哈哈,这个软件我无法安装,换了好几系统也没成功,如果是太老了,很多库文件却是。
我也懒得去解决了。
这种报错,对我这样的非计算机专业来说,简直是天书!
1
15

新的比对工具MOSAIK

功能:序列比对,类似于BWA,Bowtie
优点:全平台,甚至支持pacbio的三代测序长reads
算法:是hash index,跟其它bwt算法不太一样
作者:WP Lee - ‎2014 - ‎被引用次数:70 - ‎相关文章

Continue reading

16

NGS数据比对工具持续收集

无意中看到了这个网站,比wiki的还有全面和专业。搜集了现有还算比较出名的比对软件,并且列出来了,还做了简单评价,里面对比对工具的收集,主要是基于2012年的一个综述《Tools for mapping high-throughput sequencing data》,相信应该是有不少人都看过这篇综述的,其实生物信息初学者应该自己去文献数据库找点感兴趣的关键词的综述多看看,广泛涉猎总没有坏处的。

<img src="http://www.ebi.ac.uk/~nf/hts_mappers/mappers_timeline.jpeg" alt="Mappers Timeline" width="800">

Features Comparison

The following Table enables a comparison of mappers based on different characteristics. The table can be sorted by column (just click on the column name). The data was collected from different sources and in some cases was provided by the developers. For execution times and memory requirements we refer to the above mentioned review (supplementary data is available here).

The Data column indicates if the mapper is specifically tailored for DNA, RNA, miRNA, or bisulfite sequences.The Seq.Plat. column indicates if the mapper supports natively reads from a specific sequencing platform or not (N). The version column indicates the version of the mapper considered. Read length limits are showed in two columns: minimum read length (Min. RL) and maximum read length (Max. RL.). Unless otherwise stated the unit is base pairs. The support for mismatches and short indels is also presented including, when possible, the maximum number of allowed mismatches and indels: by default the value is presented in bases; in some cases the value is presented as a percentage of the read size; or as score, meaning that mapper uses a score function. The alignments reported column indicate the alignments reported when a read maps to multiple locations. The alignment column indicates if the reads are aligned end-to-end (Globally) or not (Locally). The Parallel column indicates if the mapper can be run in parallel and, if yes, how: using a shared-memory (SM) or/and a distributed memory (DM) computer. The QA (quality awareness) column indicates if the mapper uses read quality information during the mapping. The support for paired-end/mate-pair reads is indicated in the PE column. The Splicing column indicates, for the RNA mappers, if the detection of splice junctions is made de novo or/and through user provided libraries (Lib). The Index column indicates if the reads or/and the reference are indexed. The number of citations was obtained from Google Scholar on 13 June 2015.
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 讲的非常详细。
24

Genomemapper软件使用说明书

 我以前一直以为有了bwa跟bowtie,没什么必要用其它的alignment软件,直到我碰到了高插入删除的helicos三代测序数据,我才发现,这个古董软件genomemapper居然大有用武之地了。

一.下载并且安装该软件

这是最新版本了

Release 0.4.4 2012-10-30 source code including documentation

Wget http://1001genomes.org/data/software/genomemapper/genomemapper_0.4.4/genomemapper-0.4.4.tar.gz

这个软件安装很简单,解压进入目录,make一下即可

image001

看到make完了之后就会多了两个软件,其中一个是用来构建参考基因组索引,一个用来比对的!

二.准备数据

既然是比对软件,那么肯定是一个参考基因组,一个测序的fastq原始文件咯

当然这个软件比较奇葩,它还支持Multi-FASTA, FASTQ2 or SHORE flat file format,

三、比对命令

这里要分两步走,首先是构建参考基因组的索引,然后才是比对

/home/jmzeng/bio-soft/genomemapper-0.4.4/gmindex \

-i BRCA1.fa -x BRCA1.idx -t BRCA1.meta

首先构建索引,种子长度就用默认的12即可,然后构建完索引如下。

image002

然后进行比对即可

/home/jmzeng/bio-soft/genomemapper-0.4.4/genomemapper \

-i BRCA1.fa -q SRR258835.fastq -M 4 -G 2 -E 4 -o mapped_reads.fl -u unmapped_reads.fl

成功比对的都输出到了mapped_reads.fl -这个文件,未比对上的在unmapped_reads.fl

我有12344条序列,成功比对的只有5276条,但是如果我用精确比对的算法,只有一千五百条是可以比对的,所以用这个允许4个mismatch和2个gap的比对算法,大大提高了比对率。

然后我修改了比对参数可以达到5605,5654,5696的提升。但是没有质的飞跃,估计本身我的这种helicos测序数据错误率就太可怕了。

四,输出结果解读

image004

这个是很规则的tab键分割的文本字符,我就不解读了,大家看readme

21

SAMStat软件使用说明书

这个软件是对我们的比对结果(通常是bwa,bowtie,tophat,hisat,star)bam或者sam来进行一个可视化的总结,类似于fastqc对我们的fastq测序结果做一个可视化总结,非常好用。

一.下载并安装该软件

软件主页是http://samstat.sourceforge.net/ 里面对该软件进行非常详细的说明

包括installation和usage,我这里简单的翻译一下。

Wget http://liquidtelecom.dl.sourceforge.net/project/samstat/samstat-1.5.tar.gz

解压开看里面的readme有介绍如何安装这个软件

Unpack the tarball:

bash-3.1$ tar -zxvf samstat-XXX.tar.gz

bash-3.1$ cd samstat

bash-3.1$ ./configure

bash-3.1$ make

bash-3.1$ make check

bash-3.1$ make install

如果用root命令就可以直接用samstat啦

如果没有root权限,安装的时候稍微有点不同

./configure  --prefix=/home/jmzeng/my-bin/

make

make install

很简单的

二,数据,就是我们的bam文件啦

三,运行命令

SAMStat软件使用说明书710

四,结果

简单看看samtools flagstat 740WT1.bam  的结果

19232378 + 0 in total (QC-passed reads + QC-failed reads)

0 + 0 duplicates

18846845 + 0 mapped (98.00%:-nan%)

0 + 0 paired in sequencing

0 + 0 read1

0 + 0 read2

0 + 0 properly paired (-nan%:-nan%)

0 + 0 with itself and mate mapped

0 + 0 singletons (-nan%:-nan%)

0 + 0 with mate mapped to a different chr

0 + 0 with mate mapped to a different chr (mapQ>=5)

然后再看看我们的samstat的结果!

740WT1.bam.samstat.html

一个网页,非常丰富的内容

SAMStat软件使用说明书1168

内容太多了,我懒得解释了

见软件说明书http://davetang.org/wiki/tiki-index.php?page=SAMStat

30

Muscle进行多序列比对

软件的主页是

http://www.drive5.com/muscle/

进入主页,简单看看软件介绍,这个软件还是蛮牛的,一个人在家里自己写出来的,当然,对于普通人来说,这个软件跟clustalW没什么区别,反正都是多序列比对啦!

我们下载适合我们平台的版本即可!

Muscle进行多序列比对193

准备数据,我这里选择的是几个短小的蛋白

Muscle进行多序列比对215

 

这里有两种比对方式,都是很简单的命令

一种是先比对,再生成树文件(树的格式是Newick format, )

muscle -in mouse_J.pro -out mouse_J.pro.a

muscle -maketree -in mouse_J.pro.a -out mouse_J.phy (这里有两种构建树的方式)

另外一种是比对成aln格式的数据,然后用其它软件(phyml或者phylip)来生出树文件

muscle -in mouse_J.pro   -clwout seqs.aln

可以看到比对的效果还是蛮好的,是aln格式的比对文件,这个格式非常常用

Muscle进行多序列比对505

或者输出phy格式的比对文件,

muscle -in mouse_J.pro  -physout seqs.phy

Muscle进行多序列比对685

可以被phyml等软件识别,然后来构建进化树,见  http://www.bio-info-trainee.com/?p=626

21

Hg19基因组的分析

下载地址我就不贴了,随便谷歌一下即可!

Genome Reference Consortium Human  ---》  GRCh3

Feb. 2009 (hg19, GRCh37)这个是重点

Mar 2006 assembly = hg18 = NCBI36.

May 2004 assembly = hg17 = NCBI35.

July 2003 assembly = hg16 = NCBI34

以前的老版本就不用看啦,现在其实都已经有hg38出来啦,GRCh38 (NCBI) and hg38(UCSC)

参考:http://age.wang.blog.163.com/blog/static/119252448201092284725460/

http://www.ncbi.nlm.nih.gov/projects/genome/assembly/grc/human/

Hg19基因组的分析570

人的hg19基因组是3G的大小,因为一个英文字符是一个字节,所以也是30亿bp的碱基。

包括22条常染色体和X,Y性染色体及M线粒体染色体。

Hg19基因组的分析643

查看该文件可以看到,里面有很多的N,这是基因组里面未知的序列,用N占位,但是觉得部分都是A.T.C.G这样的字符,大小写都有,分别代表不同的意思。

然后我用linux的命令统计了一下里面这个文件的行数,

perl -lne 'END { print $. }'  hg19.fa

awk 'END { print NR }'  hg19.fa

wc -l hg19.fa

Hg19基因组的分析834

然后我写了一个脚本统计每条染色体的长度,42秒钟完成任务!

Hg19基因组的分析1125

看来这个服务器的性能还是蛮强大的,读取文件非常快!

[perl]

while(<>){

        chomp;

        if  (/>/){

if  (exists $hash_chr{$key} ){

$len = length $hash_chr{$key};

print "$key   =>   $len\n";

}

undef %hash_chr;

$key=$_;

}

else {

$hash_chr{$key}.=$_;

}

}

[/perl]

 

然后我用seed统计了一下hg19的词频(我不知道生物信息学里面的专业描述词语是什么)

Hg19基因组的分析1171

我的程序耗费了42分钟才跑完,感觉我写的程序应该是没有问题的,让我吃惊的是总共竟然只有105万条独特的10bp短序列。然后我算了一下4的10次方,(⊙o⊙)…悲剧,原来只有1048576,之所以出现这种情况,是因为里面有N这个字符串,不仅仅是A.T.C.G四个字符。我用grep -v N seed10.txt |wc -l命令再次统计了一下,发现居然就是1048576,也就是说,任意A.T.C.G四个字符组成的10bp字符串短序列在人的基因组里面都可以找到!!!

Hg19基因组的分析1407

然后我测试了一下,还是真是这样的,真是一个蛮有意思的现象。虽然我无法解释为什么,但是根据这个结果我们可以得知连续的A或者T在人类基因组里面高频出现,而连续的G或者C却很少!

如果我们储存这个10bp字符串的同时,也储存着它们在基因组的位置,那么就可以根据这个seed来进行比对,这就是blast的原理之一!

 

 

20

bowtie简单使用

首先进入bowtie的主页,千万要谷歌!!!

http://bowtie-bio.sourceforge.net/bowtie2/index.shtml

image001

主页里面有下载链接,也有索引文件,当然索引文件只有人类等模式生物

下载之后是个压缩包,解压即可使用

image003

可以看到绿色的就是命令,可以添加到环境变量使用,也可以直接用全路径使用它!!!

然后example文件夹里面有所有的测试文件。

二、准备数据

我们就用软件自带的测试数据

image005

三、运行命令

分为两步,首先索引,然后比对!!!

  • 索引,bowtie2-build your-fastq-file.fa your-index-name

image008

然后你的目录就产生了六个索引文件,我给索引取名是tmp,你们可以随便取名字

  • 然后比对,分两种,一是单端测序数据,二是双端数据

重点参数的-x 和 –S ,单端是-U 双端是-1 -2

Bowtie –x tmp –U reads.fa –S hahahhha.sam

Bowtie –x tmp -1 reads1.fa -2 reads2.fa –S hahahha.sam

四:输出文件解读

就是输出了sam文件咯,这个就看我的sam文件格式讲解哈

 

 

14

转录组比对软件tophat的使用

转录组比对软件tophat的使用

为什么要用这个软件?:因为转录组reads比对到基因组reads用bwa和bowtie的效果都不够好,所以我们选择tophat

它做了什么?:tophat把测序的转录组的原始reads比对到了参考基因组上面,并且输出了bam(二进制的sam)文件比对结果给我们。(fastq--->bam)

一:下载安装该软件

其实一般的生信服务器自然会有高手给安装好了,你只需调用即可,这里我给大家演示一下如何安装。

wget   http://ccb.jhu.edu/software/tophat/downloads/tophat-2.0.13.Linux_x86_64.tar.gz

Continue reading