[转录组学]如何执行一步式 fastq 差分分析，kallisto 来帮忙 (I)

传统转录组的分析想必大家已经非常熟悉，无非是质控 -> 比对 -> 组装 ->定量 -> 差异分析 -> 差异表达基因功能富集分析这个套路。目前公司用的流程主要也有两个门派，以“传统派”的trimmomatic/cutadapter/华大fastqc -> tophat -> cufflink -> cuffmerge -> cuffquant（较差）/RSEM -> DESeq2/edgeR -> GO/KEGG，已经较新的 fastp -> hisat2 -> stringtie -> (ballgown,较差)/ (prepDE.py+edgeR,DEseq2) -> GO/KEGG。但是对于很多样品来说我们不关注novel gene，isoform，也不关注可变剪切，是否有一套快速简单的方法一步到位的到差异分析，快速定位关键基因呢？

目的

答案当然是有的，上面讲的方法都是基于传统方法进行的分析思路。但是实际上，我们很多时候往往只是为了找到实验组和对照组的差异表达基因，以及这些基因和哪些通路有关。这个时候我们完全没有必要大费周章的对数据进行这一长串的分析。

流程

这里介绍一套轻量化的分析思路：fastp -> kallisto -> edgeR -> clusterProfiler
是的，你没有看错，所有的分析只需要这四个软件。

kallisto介绍

我们这里介绍一下这个kallisto，他是整个分析的核心软件。它是一个在2016年发表于Nature biotechnology《Near-optimal probabilistic RNA-seq quantification》的软件。
该算法认为，将fq文件的reads比对到参考基因组的具体位置，然后根据位置信息进行计数数没有实质性意义的，而该算法采用的是伪比对的方式：即直接将fq文件的reads比对参考转录组上并且直接计数。

实战

1. 准备工作

1.1 下载参考基因组

由于这里选用的是人。kallisto的使用需要我们下载物种的mRNA序列，我们可以去到ensembl去下载wget http://ftp.ensembl.org/pub/current_fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz
解压文件gunzip Homo_sapiens.GRCh38.cdna.all.fa.gz

1.2 安装软件

这里默认你已经安装好了R包和conda
conda install kallisto fastp -c bioconda -c conda-forge

1.3 下载数据

这里使用GEO的一组m6A的Input数据为例子进行分析为例，可以到https://www.ncbi.nlm.nih.gov/sra?linkname=bioproject_sra_all&from_uid=752384
进行下载。并自行搜索sra转fastq的方法，这里推荐一下多线程pfastq-dump工具。
用这个数据的原因也是埋一个坑吧，有机会我们讲一下MeRIP-seq以及和本次的RNA-seq进行关联分析。这里用到的样本是NS1-3Input和HS1-3Input，并把我们命名HS1_Input_1.fastq.gz，即样本名编号_Input_测序方向.fastq.gz

1.4 建立索引

进入到我们刚下载好的文件目录，使用命令
kallisto index -i hg38.idx Homo_sapiens.GRCh38.cdna.all.fa

2. 开始分析

2.1 准备配置信息

新建一个名为work.sh的脚本帮助我们串联流程

#!/bin/bash
index="/Index/kallisto/hg38.idx" #这里输入你的index实际路径
rawdata="/Data/GEO/PRJNA752384" #这里输入你的原始文件路径
output="/Output" #这里输入你的输出路径
cd $output #转到你的输出路径

2.2 数据质控

接下来我们就可以开始正式的分析，首先，使用fastp去接头和去低质量的碱基

mkdir cleandata
for i in $(ls $rawdata/*_1.fastq.gz) 
do
    R2=$(echo $(basename $i) | sed 's/_1/_2/g')
    sample=$(echo $(basename $i) | sed 's/_Input_1.fastq.gz//g')
    fastp -i $i -I $rawdata/$R2 -o cleandata/${sample}_1.fq.gz -O cleandata/${sample}_2.fq.gz -l 50 -q 20 -u 50 --detect_adapter_for_pe -w 8 -j cleandata/$sample.json -h cleandata/$sample.html
done

质控完成后我们可以在cleandata目录下看到每个质控的html报告。

2.3 kallisto比对

for i in $(ls cleandata/*_1.fq.gz) 
do
    R2=$(echo $i | sed 's/_1/_2/g')
    sample=$(echo $(basename $i) | sed 's/_1.fq.gz//g')
    kallisto quant -i $index -o $output/quantity/$sample -t 8 -b 100 $i $R2
done

大约过程不到5分钟每个样本，结束后我们看一下输出的文件，其中我们最关心的就是定量文件就是abundance.tsv，文件结构如下：

这里我们可以使用tpm去绘制一些热图之类的可视化分析，而其中还有一个estCount概念，他们定义如下：

即 count*(长度/有效长度)，由于DEseq2/edgeR推荐我们使用raw count所以最好在分析前进行一次转换。

本期我们就分享到这里，下一期我们会带来后续的差异分析和GO/KEGG富集分析

想了解更多？还不赶快关注我，以及我的公众号生信咖啡吗？