RNA_seq pipline
RNA_seq pipline
PeRl
2018年3月7日
首先說明一下我做RNA-seq處理流程的文件樹格式:
- RNA-seq/
- data/
- GRCh38.gtf
- chroms/
- hg38/
- samples/
- SraAccList.txt
- sra/
- fasta/
- fastqc/
- cufflinks_result/
- tophat_result/
- HTSeq_result/
- tools/
- Trimmomatic-0.36/
- data/
1. 下載參考基因組序列信息及注釋文件GTF
參考基因組用於reads的定位和識別,需要用到的文件格式是fasta以.fa為后綴, 另外的各類文件存儲格式可以參看file format.
chromFa.tar.gz是組裝后的序列,每條染色體一個文件(我們要下載的文件),用axel下載數據
#例如我們要下載hg38,其中-n 20表示線程數
axel -n 20 http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.chromFa.tar.gz解壓后,我們需要把染色體文件整合成一個整體文件,然后用bowtie2建立索引文件
#! /usr/bin/env bash
for i in `ls /home/wang/Downloads/RNA-seq/data/chroms `
do
`cat /home/wang/Downloads/RNA-seq/data/chroms/$i >> /home/wang/Downloads/RNA-seq/data/hg38/hg38.fa`
done建立索引文件(有點耗時,約2小時左右):
bowtie2-build /home/test/RNA-seq/data/hg38/hg38.fa hg38從UCSC下載注釋文件:
2..sra數據下載及格式轉換
.sra數據的處理需要用到sra-toolkit
sudo apt-get install sra-toolkit然后用Bash腳本批量下載數據
#! /usr/bin/env bash
for i in `cat /home/wang/Downloads/RNA-seq/data/samples/SraAccList.txt`
do
echo $i
`axel -n 20 ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/ERR/ERR212/$i/$i.sra > /home/wang/Downloads/RNA-seq/data/samples/sra/$i.sra`
# `mkdir /home/wang/Downloads/RNA-seq/data/samples/fasta/$i`
# `fastq-dump -O /home/wang/Downloads/RNA-seq/data/samples/fasta/$i --split-3 /home/wang/Downloads/RNA-seq/data/samples/sra/$i.sra`
done得到樣本的.sra數據后,我們仍需要轉換為.fa數據:
--split-3: 拆分文件,如果得到的.sra文件是單末端,那么這個參數就會被忽略;如果原文件中有兩個文件,那么它就會把成對的文件按*_1.fastq,*_2.fastq這樣分開。
fastq-dump [-O <outdir path>] –split-3 <input sra path>
fastq-dump -O /home/wang/Downloads/RNA-seq/data/samples/fasta/ERR2124441 --split-3 /home/wang/Downloads/RNA-seq/data/samples/sra/ERR2124441.sra3. Fasta質量控制
FastQC最為重要的部分是：Per base sequence quality（鹼基質量）,Per base sequnence content（鹼基含量分布）,Adapter Content（接頭含量）。
sudo apt-get install fastqc
fastqc -q /home/wang/Downloads/RNA-seq/data/samples/fasta/ERR2124441/*.fastq -o /home/wang/Downloads/RNA-seq/data/samples/fastqc- 總體的質量表現
上圖中Summary的部分就是整個報告的目錄,整個報告分成若干個部分。合格會有個綠色的對勾,警告是個“!”,不合格是個紅色的叉子。
- 鹼基質量
一般認為從第二個鹼基開始,平均每個鹼基的測序質量在四分位線在30分以上,則認為測序質量非常好,也就是我們通常稱的Q30。
- 每條序列的測序質量統計
統計序列測序質量,如果均值在高分段就說明質量可以。
- 鹼基含量統計
統計每種鹼基含量,一般前面幾bp序列可能包含引物而導致質量不高,可以cut掉；后面部分應該在25%左右,像下面的例子就是GC含量不合格。
- GC含量分布
藍色的線是程序根據經驗分布給出的理論值,紅色是真實值,兩個應該比較接近才比較好
- adapter
如果有adapter序列沒有去除干凈的情況,在后續分析的時候需要先使用Trimmomatic軟件。
Trimmomatic需要區分PE(Paired End)和SE(Single End).
java -jar <path to trimmomatic.jar> PE/SE [-threads <threads] [-phred33 | -phred64] [-trimlog <logFile>] <input 1> <input 2> <paired output 1> <unpaired output 1> <paired output 2> <unpaired output 2> <step 1> …
過濾數據的步驟與命令行中過濾參數的順序有關,參數之間用空格分隔,參數和參數數值之間用:分隔,通常的過濾步驟如下：
| 參數 | 意義 | options |
|---|---|---|
| ILLUMINACLIP | 過濾reads中的Illumina測序接頭和引物序列,並決定是否去除反向互補的R1/R2中的R2 | <需要的adapter地址>:<最大錯配數>:<跨adapter序列的匹配精確度>:<兩個adapter之間序列的匹配精確度> |
| SLIDINGWINDOW | 從 reads 的 5’ 端開始,進行滑窗質量過濾,切掉鹼基質量平均值低於閾值的滑窗 | <滑動窗大小>:<質量平均值> |
| MAXINFO | 一個自動調整的過濾選項,在保證 reads 長度的情況下盡量降低測序錯誤率,最大化 reads 的使用價值 | NULL |
| LEADING | 從 reads 的開頭切除質量值低於閾值的鹼基 | <最小質量> |
| TRAILING | 從 reads 的末尾開始切除質量值低於閾值的鹼基 | <最小質量> |
| CROP | 從 reads 的末尾切掉部分鹼基使得 reads 達到指定長度 | <長度> |
| HEADCROP | 從 reads 的開頭切掉指定數量的鹼基 | <鹼基數量> |
| MINLEN | 如果經過剪切后 reads 的長度低於閾值則丟棄這條 reads | <最短長度> |
| AVGQUAL | 如果 reads 的平均鹼基質量值低於閾值則丟棄這條 reads | NULL |
| TOPHRED33 | 將 reads 的鹼基質量值體系轉為 phred-33 | NULL |
| TOPHRED64 | 將 reads 的鹼基質量值體系轉為 phred-64 | NULL |
java -jar /home/wang/Downloads/RNA-seq/tools/Trimmomatic-0.36/trimmomatic-0.36.jar PE -threads 12 -phred33 -trimlog /home/wang/Downloads/RNA-seq/tools/Trimmomatic-0.36/logfile.txt /home/wang/Downloads/RNA-seq/data/samples/fasta/ERR2124441/ERR2124441_1.fastq /home/wang/Downloads/RNA-seq/data/samples/fasta/ERR2124441/ERR2124441_2.fastq /home/wang/Downloads/RNA-seq/data/samples/fasta/ERR2124441/ERR2124441_1_paired.fq /home/wang/Downloads/RNA-seq/data/samples/fasta/ERR2124441/ERR2124441_1_unpaired.fq /home/wang/Downloads/RNA-seq/data/samples/fasta/ERR2124441/ERR2124441_2_paired.fq /home/wang/Downloads/RNA-seq/data/samples/fasta/ERR2124441/ERR2124441_2_unpaired.fq ILLUMINACLIP:/home/wang/Downloads/RNA-seq/tools/Trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:504. Tophat mapping
處理完reads后，我們需要用tophat把reads匹配到參考基因組上。
tophat -o <output path> -G <gtf path> <索引文件path> <fasta文件path>
tophat -o /home/wang/Downloads/RNA-seq/data/samples/tophat_result/ERR2124441 -G /home/wang/Downloads/RNA-seq/data/GRCh38.gtf /home/wang/Downloads/RNA-seq/data/hg38/hg38 /home/wang/Downloads/RNA-seq/data/samples/fasta/ERR2124441/*_paired.fq最后得到的結果如下，一般情況下，接下去需要用到的就是accepted_hits.bam
接着我們需要檢查一下得到的accepted_hits.bam是否已經排序:
sudo apt-get install samtools
samtools view -h /home/wang/Downloads/RNA-seq/data/samples/tophat_result/ERR2124441/accepted_hits.bam | head如果header中顯示SO:coordinate,那就說明已經排序了。
5. RNA-Seq 質控
基於原始測序數據的質控（例如FastQC）不足以保證RNA-seq 數據的可用性。
我們還需要使用RSeQC程序包(基於python)全面地評估RNA-seq 結果質量，例如序列質量、GC 偏倚、PCR 偏倚、核苷酸組成偏倚、序列深度、鏈特異性、覆蓋均一性，和基因組結構上的片段分布等，以確保后續分析的可靠性。
#安裝RSeQc
sudo apt-get install python-pip
pip install RSeQc
#如果出現compr_lzo.c:31:23: fatal error: lzo/lzo1x.h: No such file or directory,需要安裝 liblzo2-dev
sudo apt-get install liblzo2-dev簡單的用bam_stat.py看下統計信息
[bam_stat.py的路徑] -i [.bam的路徑] > [輸出路徑]
#找到python安裝包的文件夾地址
.local/bin/bam_stat.py -i /home/wang/Downloads/RNA-seq/data/samples/tophat_result/ERR2124441/accepted_hits.bam > /home/wang/Downloads/RNA-seq/data/samples/tophat_result/ERR2124441/output_bam_stat.txt
6. cufflinks得到fpkm值
因為前面用的GRch38.gtf用的是Refseq ID,最后還需要ID轉換，這個就不贅述了。
#sudo apt-get install cufflinks
cufflinks -p 4 -G /home/wang/Downloads/RNA-seq/data/GRCh38.gtf -o /home/wang/Downloads/RNA-seq/data/samples/cufflinks_result/ERR2124441 /home/wang/Downloads/RNA-seq/data/samples/tophat_result/ERR2124441/accepted_hits.bam7. HTSeq得到count數及DEseq2處理
處理count數用到的是python工具HTSeq
pip install HTSeqhtseq-count使用時需用samtools對accepted_hist.bam文件以read name排序（-n）才能使用
samtools sort -n /home/wang/Downloads/RNA-seq/data/samples/tophat_result/ERR2124441/accepted_hits.bam -o /home/wang/Downloads/RNA-seq/data/samples/tophat_result/ERR2124441/accepted_hits_sort.bam 排序后用htseq-count計數,輸入文件為排序后的tophat比對結果accepted_hits_sort.bam
根據不同需要可以修改參數，參看官方說明,其中-i gene_id描述的是使用的:
/home/wang/.local/bin/htseq-count -f bam -s no -a 10 -t exon -i gene_id -m union /home/wang/Downloads/RNA-seq/data/samples/tophat_result/ERR2124441/accepted_hits_sort.bam /home/wang/Downloads/RNA-seq/data/samples/cufflinks_result/ERR2124441/transcripts.gtf > /home/wang/Downloads/RNA-seq/data/samples/HTSeq_result/ERR2124441/SRR1206035.htseq_count
接下去的處理就是用DESeq2處理count數，以前我有做過總結，這方面就不再加以說明了，就是醬紫~\(≧▽≦)/~