傳統(tǒng)轉(zhuǎn)錄組的分析想必大家已經(jīng)非常熟悉，無非是
質(zhì)控 -> 比對 -> 組裝 ->定量 -> 差異分析 -> 差異表達基因功能富集分析這個套路。目前公司用的流程主要也有兩個門派，以“傳統(tǒng)派”的trimmomatic/cutadapter/華大fastqc -> tophat -> cufflink -> cuffmerge -> cuffquant（較差）/RSEM -> DESeq2/edgeR -> GO/KEGG，已經(jīng)較新的 fastp -> hisat2 -> stringtie -> (ballgown,較差)/ (prepDE.py+edgeR,DEseq2) -> GO/KEGG。但是對于很多樣品來說我們不關(guān)注novel gene，isoform，也不關(guān)注可變剪切，是否有一套快速簡單的方法一步到位的到差異分析，快速定位關(guān)鍵基因呢？
目的
答案當然是有的，上面講的方法都是基于傳統(tǒng)方法進行的分析思路。但是實際上，我們很多時候往往只是為了找到實驗組和對照組的差異表達基因，以及這些基因和哪些通路有關(guān)。這個時候我們完全沒有必要大費周章的對數(shù)據(jù)進行這一長串的分析。
流程
這里介紹一套輕量化的分析思路：fastp -> kallisto -> edgeR -> clusterProfiler
是的，你沒有看錯，所有的分析只需要這四個軟件。
kallisto介紹
我們這里介紹一下這個kallisto，他是整個分析的核心軟件。它是一個在2016年發(fā)表于Nature biotechnology《Near-optimal probabilistic RNA-seq quantification》的軟件。
該算法認為，將fq文件的reads比對到參考基因組的具體位置，然后根據(jù)位置信息進行計數(shù)數(shù)沒有實質(zhì)性意義的，而該算法采用的是偽比對的方式：即直接將fq文件的reads比對參考轉(zhuǎn)錄組上并且直接計數(shù)。
實戰(zhàn)
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 安裝軟件
這里默認你已經(jīng)安裝好了R包和conda
conda install kallisto fastp -c bioconda -c conda-forge
1.3 下載數(shù)據(jù)
這里使用GEO的一組m6A的Input數(shù)據(jù)為例子進行分析為例，可以到https://www.ncbi.nlm.nih.gov/sra?linkname=bioproject_sra_all&from_uid=752384
進行下載。并自行搜索sra轉(zhuǎn)fastq的方法，這里推薦一下多線程pfastq-dump工具。
用這個數(shù)據(jù)的原因也是埋一個坑吧，有機會我們講一下MeRIP-seq以及和本次的RNA-seq進行關(guān)聯(lián)分析。這里用到的樣本是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的腳本幫助我們串聯(lián)流程
#!/bin/bash
index="/Index/kallisto/hg38.idx" #這里輸入你的index實際路徑
rawdata="/Data/GEO/PRJNA752384" #這里輸入你的原始文件路徑
output="/Output" #這里輸入你的輸出路徑
cd $output #轉(zhuǎn)到你的輸出路徑
2.2 數(shù)據(jù)質(zhì)控
接下來我們就可以開始正式的分析，首先，使用fastp去接頭和去低質(zhì)量的堿基
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
質(zhì)控完成后我們可以在cleandata目錄下看到每個質(zhì)控的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分鐘每個樣本，結(jié)束后我們看一下輸出的文件，其中我們最關(guān)心的就是定量文件就是abundance.tsv，文件結(jié)構(gòu)如下：
這里我們可以使用tpm去繪制一些熱圖之類的可視化分析，而其中還有一個
estCount概念，他們定義如下：
即 count*(長度/有效長度)，由于DEseq2/edgeR推薦我們使用raw count所以最好在分析前進行一次轉(zhuǎn)換。
本期我們就分享到這里，下一期我們會帶來后續(xù)的差異分析和GO/KEGG富集分析
想了解更多？還不趕快關(guān)注我，以及我的公眾號生信咖啡嗎？
