一考赛、測序數(shù)據(jù)的格式轉(zhuǎn)換
sra文件下載好后,使用fastq-dump轉(zhuǎn)換數(shù)據(jù)格式:
fastq-dump --split-files SRR6232298.sra
fastq-dump --gzip --split-files SRR6232298.sra
#轉(zhuǎn)換格式的同時解壓為gz文件莉测,節(jié)省空間
當(dāng)有多個sra文件時颜骤,可通過腳本進(jìn)行批量解壓:
vi fastq-dump.sh? ? #創(chuàng)建一個腳本文件,內(nèi)容如下:
#捣卤!/bin/bash
for i in ~ncbi/public/sra/sra*
do
echo $i
fastq-dump --gzip --split-files $i
done
echo OK
二忍抽、測序數(shù)據(jù)的質(zhì)控
FastQC---測序數(shù)據(jù)質(zhì)控的軟件?
是一個java軟件,下載后可以直接使用(免安裝編譯)董朝,但是需要自行配置好java環(huán)境
首先我們配置java環(huán)境(已下好java文件梯找,為下述的jdk-8u172-linux-x64.tar.gz):
sudo mkdir /usr/java
sudo tar -zvxf /home/noodles/Biosofts/jdk-8u172-linux-x64.tar.gz -C /usr/java/
sudo cd /usr/java
sudo ln -s jdk1.8.0_172 latest
sudo ln -s /usr/java/latest default
sudo vi /etc/profile
末尾加上如下三行:?
export JAVA_HOME=/usr/java/latest
export PATH=$JAVA_HOME/bin:$JAVA_HOME/jre/bin:$PATH
export CLASSPATH=.:$JAVA_HOME/lib/dt.jar:$JAVA_HOME/lib/tools.jar
source /etc/profile
java -version
配置好java環(huán)境后,接著可以開始下載安裝FastQC軟件了
1. wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.7.zip?
2. unzip /home/noodles/Biosofts/fastQC/fastqc_v0.11.7.zip -d ~/Biosofts/fastQC
3. ~/Biosofts/fastQC/FastQC/fastqc -h
當(dāng)我運(yùn)行這步后益涧,系統(tǒng)提示fastqc還未安裝,通過apt install fastqc來安裝驯鳖,當(dāng)進(jìn)行安裝后報(bào)錯:
呈上:
4. 將fastqc加入環(huán)境變量:
echo 'export PATH=~/Biosofts/fastQC/FastQC:$PATH' >>~/.bashrc
source ~/.bashrc
fastqc -h
至此御蒲,已裝好了fastQC軟件,fastQC的使用方法如下:
fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN
各參數(shù)的含義:
-o --outdir 生成的報(bào)告文件的儲存路徑诊赊,生成的報(bào)告的文件 名是根據(jù)輸入來定的
--extract 生成的報(bào)告默認(rèn)會打包成1個壓縮文件厚满,使用這個 參數(shù)是讓程序不打包?
-t --threads 選擇程序運(yùn)行的線程數(shù)?
-c --contaminants 污染序列選項(xiàng),輸入的是一個文件碧磅,格式 是Name [Tab] Sequence碘箍,里面是可能的污染序列
-a --adapters 也是輸入一個文件,文件的格式Name [Tab] Sequence鲸郊,儲存的是測序的adpater序列信息丰榴;默認(rèn)已有一些?
-q --quiet 安靜運(yùn)行模式,一般不選這個選項(xiàng)的時候秆撮,程序 會實(shí)時報(bào)告運(yùn)行的狀況四濒。?
接下來以用fastqc處理之前下載的seq為例,當(dāng)我處理SRR6232298_1.fastq.gz時報(bào)錯:
經(jīng)過檢查后原來是fastqc的權(quán)限問題,增加權(quán)限即可:
三复隆、測序數(shù)據(jù)過濾
數(shù)據(jù)過濾軟件用來切除接頭序列和低質(zhì)量堿基,目前也已有很多工具:Trimmomatic姆涩、seqtk挽拂、 cutadapt、 bbduk(BBmap). 下面以Trimmomatic為例介紹:
Trimmomatic 是一個廣受歡迎的 Illumina 平臺數(shù)據(jù)過濾工具骨饿。 Trimmomatic 支持多線程亏栈,處理數(shù)據(jù)速度快,主要用來去除 Illumina 平臺的 Fastq 序列中的接頭宏赘,并根據(jù)堿基質(zhì)量值對 Fastq 進(jìn)行修剪绒北。軟件有兩種過濾模式,分別對應(yīng) SE 和 PE 測序數(shù)據(jù)察署,同時支持 gzip 和 bzip2 壓縮文件闷游。另外也支持 phred-33 和 phred-64 格式互相轉(zhuǎn)化,不過現(xiàn)在 絕大部分 Illumina 平臺的產(chǎn)出數(shù)據(jù)也都轉(zhuǎn)為使用 phred-33 格式了 贴汪。
下載:
wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.38.zip
安裝:?
unzip Trimmomatic-0.38.zip -d ~/Biosofts/trimmomatic/Trimmomatic038/?
運(yùn)行:?
java -jar ~/Biosofts/trimmomatic/Trimmomatic038/Trimmomatic-0.38/trimmomatic-0.38.jar?
添加環(huán)境變量:
echo 'export PATH=~/Biosofts/trimmomatic/Trimmomatic038/Trimmomatic-0.38/trimmomatic-0.38.jar:$PATH' >>~/.bashrc
source ~/.bashrc
然后以之前下載的SRR6232298_1.fastq.gz數(shù)據(jù)為例進(jìn)行操作:
java -jar ~/Biosofts/trimmomatic/Trimmomatic038/Trimmomatic-0.38/trimmomatic-0.38.jar PE -phred33?SRR6232298_1.fastq.gz?./trim_out/output_forward_paired.fq.gz ./trim_out/output_forward_unpaired.fq.gz ./trim_out/output_reverse_paired.fq.gz ./trim_out/output_reverse_unpaired.fq.gz?ILLUMINACLIP:/home/noodles/Biosofts/trimmomatic/Trimmomatic038/Trimmomatic-0.38/adapters/TruSeq2-PE.fa:2:30:10?SLIDINGWINDOW:5:20 LEADING:20 TRAILING:20 MINLEN:75
(上述命令中脐往,output_forward_paired.fq.gz、output_forward_unpaired.fq.gz扳埂、output_reverse_paired.fq.gz业簿、output_reverse_unpaired.fq.gz四個文件由自己命名,軟件自行生成阳懂。)結(jié)果報(bào)錯:
經(jīng)仔細(xì)檢查后岩调,發(fā)現(xiàn)問題出在下面的點(diǎn):
trimmomatic有兩種模式即SE模式和PE模式巷燥,分別對應(yīng)單末端測序模式和雙末端測序模式,在 SE 模式下号枕,只有一個輸入文件和一個過濾之后的輸出文件矾湃;在 PE 模式下,有兩個輸入文件堕澄,正向測序序列和反向測序序列邀跃,但是過濾之后輸出文件有四個,過濾之后雙端序列都保留的就是 paired蛙紫,反之如果其中一端序列過濾之后被丟棄了另一端序列保留下來了就是 unpaired 拍屑。上述的問題就在于選用的是PE模式,但只輸入了一個文件坑傅,所以報(bào)錯僵驰。(SRR6232298.sra為雙末端測序的數(shù)據(jù),SRR6232298_1.fastq.gz和SRR6232298_2.fastq.gz為使用fastq-dump對SRR6232298.sra處理后所得,分別為正向測序序列和反向測序序列)因此當(dāng)我將SRR6232298_1.fastq.gz和SRR6232298_2.fastq.gz都輸入時蒜茴,軟件可正常運(yùn)行:
java -jar ~/Biosofts/trimmomatic/Trimmomatic038/Trimmomatic-0.38/trimmomatic-0.38.jar PE -phred33 SRR6232298_1.fastq.gz SRR6232298_2.fastq.gz ./trim_out/output_forward_paired.fq.gz ./trim_out/output_forward_unpaired.fq.gz ./trim_out/output_reverse_paired.fq.gz ./trim_out/output_reverse_unpaired.fq.gz ILLUMINACLIP:/home/noodles/Biosofts/trimmomatic/Trimmomatic038/Trimmomatic-0.38/adapters/TruSeq2-PE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:20 TRAILING:20 MINLEN:75?
Trimmomatic 過濾數(shù)據(jù)的步驟與命令行中過濾參數(shù)的順序有關(guān)星爪,通常的過濾步驟如下:
1.?ILLUMINACLIP: 過濾 reads 中的 Illumina 測序接頭和引物序列,并決定是否去除反向互補(bǔ)的 R1/R2 中的 R2粉私。
2.?SLIDINGWINDOW: 從 reads 的 5' 端開始顽腾,進(jìn)行滑窗質(zhì)量過濾,切掉堿基質(zhì)量平均值低于閾值的滑窗诺核。
3.?MAXINFO: 一個自動調(diào)整的過濾選項(xiàng)抄肖,在保證 reads 長度的情況下盡量降低測序錯誤率,最大化 reads 的使用價(jià)值窖杀。
4.?LEADING: 從 reads 的開頭切除質(zhì)量值低于閾值的堿基漓摩。
5.?TRAILING: 從 reads 的末尾開始切除質(zhì)量值低于閾值的堿基。
6.?CROP: 從 reads 的末尾切掉部分堿基使得 reads 達(dá)到指定長度入客。
7.?HEADCROP: 從 reads 的開頭切掉指定數(shù)量的堿基管毙。
8.?MINLEN: 如果經(jīng)過剪切后 reads 的長度低于閾值則丟棄這條 reads。
9. AVGQUAL: 如果 reads 的平均堿基質(zhì)量值低于閾值則丟棄這條 reads桌硫。
10.?TOPHRED33: 將 reads 的堿基質(zhì)量值體系轉(zhuǎn)為 phred-33夭咬。
11.?TOPHRED64: 將 reads 的堿基質(zhì)量值體系轉(zhuǎn)為 phred-64。
參考:NGS 數(shù)據(jù)過濾之 Trimmomatic 詳細(xì)說明
下面兩張圖分別為未經(jīng)Trimmomatic過濾和經(jīng)過Trimmomatic過濾的SRR6232298_1.fastq.gz鞍泉,經(jīng)過fastQC質(zhì)控后的結(jié)果: