Kallisto使用说明.docx

Kallisto使用说明2017-03-19 Ren NanFang HospitalKallisto是一款快速比对的软件，由加州大学伯克利分校的Bray等人于2015年开发出来。它的优势是快速、耗费内存小，可以在普通台式机上在几分钟之内完成人或其他物种的转录组二代测序比对任务。其最大的特点是不需要完整的参考基因组，但是也需要对某物种的全转录本序列建立索引，然后再进行所谓的假比对（pseudoalignment）。该算法的基础是转录组德布鲁因图（transcriptome de Bruijn graph， T-DBG）Kallisto可以从github上获得二进制文件，https:/pachterlab.github.io/kallisto/，目前的版本是v0.43.0，2016年6月2日更新。在Linux系统上对文件进行解压缩，进入相应文件夹kallisto_linux-v*后可以看到可执行文件kallisto，把该文件的路径写入到PATH中，即可在其它路径下运行。如下图所示，在用户目录下键入kallisto则出现软件的版本和简单的使用说明。在github网页上有Kallisto的使用手册。从上图可以看到Kallisto主要有6个命令，分别是index，quant，pseudo，h5dump，version，cite。其中最常使用的是前2个，index建立转录组索引和quant进行转录本水平的表达定量。pseudo为单细胞转录组测序使用，h5dump为转换文件格式，把hdf5格式的文件转换成纯文本tsv格式的文件。另外2个version则是查阅软件版本号，cite则是显示出引用文献。kallisto 0.43.0Usage: kallisto arguments Where can be one of: index Builds a kallisto index quant Runs the quantification algorithm pseudo Runs the pseudoalignment step h5dump Converts HDF5-formatted results to plaintext version Prints version information cite Prints citation informationRunning kallisto without arguments prints usage information for 所需材料： Linux电脑（我是用的版本是Ubuntu 16.04 Sever LTS x64）、Kallisto、测序文件Fastaq格式或者Fastaq压缩文件、相应物种的全转录组cdna序列，可以从UCSC或者Ensembl上下载获得，我是用的是Ensembl的cdna序列文件。假设以上材料已经准备好，且Kallisto安装完可以运行，那么下面要做的第一件事就是使用该物种的全转录组cdna序列文件建立索引，使用kallisto index命令，如下：kallisto index PATH/Mus_musculus.GRCm38.cdna.all.fa.gz -i ensemblgrc38.mm.87输入相应cdna序列的压缩文件，-i后跟建立的索引文件的名字。几分钟后可以看到文件夹中有了ensemblgrc38.mm.87的索引文件。第二，对测序数据进行定量分析使用kallisto quant命令，这里我使用的是小鼠10.5天胚胎的RNA-seq双端测序数据（*.R1.clean.fastq.gz & *.R2.clean.fastq.gz），如下：kallisto quant -i PATH/ensemblgrc38.mm.87 -o ./ -t 24 PATH/10-5_L2_I367.R1.clean.fastq.gz PATH/10-5_L2_I367.R2.clean.fastq.gz参数-i后面跟刚才建立好的索引名字，-o后面跟输出文件的指定目录；其它选项如-t后面跟线程数量，最后是FASTAQ格式的测序文件，如果是双端测序（paired-end）文件两个文件中间使用空格，程序会自动识别两个文件为左右两端的测序文件。如果使用的是单端测序（single-end）文件，要使用-single参数指定；并且一定要使用-l参数，后面跟片段长度，估计的平均片段长度；-s参数后面跟估计的片段长度的标准差。-l和-s参数最好使用类似Agilent Bioanalyzer软件去确定。kallisto quant -i index -o output -single -l 200 -s 20 file1.fastq.gz file2.fastq.gz file3.fastq.gz约35分钟后，在相应文件夹就会得到3个输出文件：abundances.h5、abundances.tsv和run_info.json。abundances.h5： HDF5二进制格式的文件，包含了运行日志信息、表达丰度估计值、bootstrap估计和转录本长度信息。该文件可以直接用sleuth读取处理，也可以使用kallisto h5dump命令将其转变为纯文本的tsv格式文件；abundances.tsv： 包含有表头的纯本文tsv格式文件，表头是：target_id, length, eff_length, est_counts, tpm；run_info.json： 一个json格式的日志文件。下游分析：得到转录本表达定量文件后，一般来说下游分析是要寻找差异表达基因。与kallisto搭档的下游分析软件是sleuth，如果我们想要使用R的Bioconductor中的分析包比如DESeq2或者edgeR的话，那么就需要使用tximport软件包，将kallisto得到的文件读入R中进行处理后，基因表达counts值可以被DESeq2或者edgeR处理。下面展示如何把kallisto生成的文件转换成被R包读取的文件，说明摘自/packages/release/bioc/vignettes/tximport/inst/doc/tximport.html 网页。首先安装常用软件并加载：source(/biocLite.R)biocLite(edgeR)biocLite(DESeq2)biocLite(tximport)biocLite(readr)library(tximport)library(edgeR)library(DESeq2)加载这些包之后，我们需要该物种的个转录本名称和基因名称的对照数据，如下面所示： tx_id symbol1 ENSMUST00000000001 Gnai32 ENSMUST00000000003 Pbsn3 ENSMUST00000114041 Pbsn4 ENSMUST00000000028 Cdc455 ENSMUST00000115585 Cdc456 ENSMUST00000096990 Cdc457 ENSMUST00000152754 H198 ENSMUST00000140716 H199 ENSMUST00000132294 H1910 ENSMUST00000149974 H19这个对照表可以从EnsDB包中获取，添加EnsDb.Mmusculus.v79，ensembldb包，如果是人类基因组的话那就下载EnsDb.Hsapiens.v86包。加载包之后，通过以下命令得到基因和转录本名称对照表。library(EnsDb.Mmusculus.v79)library(ensembldb)tx.mms - EnsDb.Mmusculus.v79tx.mms - transcripts(mmsdb, return.type = DataFrame)tx.mms - tx.mms, c(tx_id, gene_id)gene.mms - genes(mmsdb, return.type = DataFrame)gene.mms - gene.mms, c(gene_id, symbol)tx2genes.mms - merge(tx.mms, gene.mms, by = gene_id)tx2genes - tx2genes.mms, 2:3head(tx2genes) tx_id symbol1 ENSMUST00000000001 Gnai32 ENSMUST00000000003 Pbsn3 ENSMUST00000114041 Pbsn4 ENSMUST00000000028 Cdc455 ENSMUST00000115585 Cdc456 ENSMUST00000096990 Cdc45这里我们有5个小鼠胚胎样本的RNA-seq数据，经过以上kallisto定量分析后，我们得到5个abundance.csv 使用tximport命令将这5个文件读入R，在type参数里指明数据是kallisto生成的文件，type参数里还可以指定其它软件生成的文件如salmon、sailfish、rsem。txi - tximport(files = files, type = kallisto, tx2gene = tx2gene, reader = read_tsv)names(txi)1 abundancecountslengthcountsFromAbundancehead(txi$counts) 7.5d 10.5d 13.5d 16.5d 18.5d 3512.0930 2904.7000 2807.9760 9215.1200 8124.41000610007P14Rik1103.4468 655.5391 1130.1470 1212.0528 637.99970610009B22Rik226.0001 168.0000 141.0000 192.0000 302.00000610009L18Rik11.0000 19.0000 27.0000 36.0000 11.00000610009O20Rik1339.0000 772.0000 1460.0000 2202.0000 1007.00000610010F05Rik402.9997 291.9996 445.9999 532.9993 338.9997我们也可以使用countsFromAbundance参数中的选项，对估计的counts值进行尺度转变 (scaled)，scaledTPM或者lengthScaledTPM，转变后counts值就不再与转录本长度有关。使用这两个方法之后，我们就不用再为下游分析软件提供库大小或者转录本长度的信息了。具体使用方法参见tximport帮助文档。我们同样可以使用txOut参数，指定是否进行转录本水平的表达分析，而非基因水平的表达分析，默认值是FALSE，即对同一基因的所有转录本的counts进行合并来代表该基因的counts值，而TRUE则给出原始转录本水平的counts值估计。将Kallisto定量软件的转录本丰度估计值导入edgeR、DESeq或者limma-voom等分析软件。1. 下面是edgeR的数据输入方法：cts - txi$countsnormMat - txi$lengthnormMat - normMat/exp(rowMeans(log(normMat)library(edgeR)o - log(calcNormFactors(cts/normMat) + log(colSums(cts/normMat)y - DGEList(cts)y$offset - t(t(log(normMat) + o)# y is now ready for estimate dispersion functions see edgeR Users Guide2. 以上得到的y就可以用于估计散度系数并进行后续的分析了。下面看一下使用DESeq2分析时需要什么样的数据输入。使用者最好是先建立一个样本表格sampletale，sampletable的行名字要与txi$counts中的列名字相对应，即样本的名称。这里我们创建一个sampletable它只含有各样本的处理条件，我们这里使用的是简单时间序列，那么对于每一个小鼠样本的时间点就是每个样本的条件。sampletable - data.frame(condition = factor(c(7.5d, 10.5d, 13.5d, 16.5d, 18.5d)rownames(sampletable) - colnames(txi$counts)sampletable# condition# 7.5d 7.5d# 10.5d 10.5d# 13.5d 13.5d# 16.5d 16.5d# 18.5d 18.5ddds