測序數(shù)據(jù)的格式轉(zhuǎn)換與質(zhì)控

一考赛、測序數(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)錯:

當(dāng)按照提示apt-get update后再進(jìn)行安裝即可闲询。此命令的兩個作用:1、apt-get?update是同步 /etc/apt/sources.list 和 /etc/apt/sources.list.d 中列出的源的索引浅辙,這樣才能獲取到最新的軟件包扭弧。2、apt-get?update只是更新了apt的資源列表记舆,沒有真正的對系統(tǒng)執(zhí)行更新鸽捻。如果需要,要使用apt-get?upgrade來更新泽腮。

呈上:

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)限即可:

fastqc可一條命令對多個序列進(jìn)行指控峻黍,也可以通過腳本進(jì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)錯:

但是確認(rèn)多次adapter文件路徑后梅尤,仍然報(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é)果:

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末,一起剝皮案震驚了整個濱河市肮帐,隨后出現(xiàn)的幾起案子咖驮,更是在濱河造成了極大的恐慌,老刑警劉巖训枢,帶你破解...
    沈念sama閱讀 216,324評論 6 498
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件托修,死亡現(xiàn)場離奇詭異,居然都是意外死亡恒界,警方通過查閱死者的電腦和手機(jī)睦刃,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 92,356評論 3 392
  • 文/潘曉璐 我一進(jìn)店門,熙熙樓的掌柜王于貴愁眉苦臉地迎上來十酣,“玉大人涩拙,你說我怎么就攤上這事∷什桑” “怎么了兴泥?”我有些...
    開封第一講書人閱讀 162,328評論 0 353
  • 文/不壞的土叔 我叫張陵,是天一觀的道長虾宇。 經(jīng)常有香客問我搓彻,道長,這世上最難降的妖魔是什么? 我笑而不...
    開封第一講書人閱讀 58,147評論 1 292
  • 正文 為了忘掉前任旭贬,我火速辦了婚禮怔接,結(jié)果婚禮上,老公的妹妹穿的比我還像新娘稀轨。我一直安慰自己扼脐,他們只是感情好,可當(dāng)我...
    茶點(diǎn)故事閱讀 67,160評論 6 388
  • 文/花漫 我一把揭開白布靶端。 她就那樣靜靜地躺著谎势,像睡著了一般。 火紅的嫁衣襯著肌膚如雪杨名。 梳的紋絲不亂的頭發(fā)上脏榆,一...
    開封第一講書人閱讀 51,115評論 1 296
  • 那天,我揣著相機(jī)與錄音台谍,去河邊找鬼须喂。 笑死,一個胖子當(dāng)著我的面吹牛趁蕊,可吹牛的內(nèi)容都是我干的坞生。 我是一名探鬼主播,決...
    沈念sama閱讀 40,025評論 3 417
  • 文/蒼蘭香墨 我猛地睜開眼掷伙,長吁一口氣:“原來是場噩夢啊……” “哼是己!你這毒婦竟也來了?” 一聲冷哼從身側(cè)響起任柜,我...
    開封第一講書人閱讀 38,867評論 0 274
  • 序言:老撾萬榮一對情侶失蹤卒废,失蹤者是張志新(化名)和其女友劉穎,沒想到半個月后宙地,有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體摔认,經(jīng)...
    沈念sama閱讀 45,307評論 1 310
  • 正文 獨(dú)居荒郊野嶺守林人離奇死亡,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 37,528評論 2 332
  • 正文 我和宋清朗相戀三年宅粥,在試婚紗的時候發(fā)現(xiàn)自己被綠了参袱。 大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片。...
    茶點(diǎn)故事閱讀 39,688評論 1 348
  • 序言:一個原本活蹦亂跳的男人離奇死亡秽梅,死狀恐怖抹蚀,靈堂內(nèi)的尸體忽然破棺而出,到底是詐尸還是另有隱情企垦,我是刑警寧澤况鸣,帶...
    沈念sama閱讀 35,409評論 5 343
  • 正文 年R本政府宣布,位于F島的核電站竹观,受9級特大地震影響镐捧,放射性物質(zhì)發(fā)生泄漏潜索。R本人自食惡果不足惜,卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 41,001評論 3 325
  • 文/蒙蒙 一懂酱、第九天 我趴在偏房一處隱蔽的房頂上張望竹习。 院中可真熱鬧,春花似錦列牺、人聲如沸整陌。這莊子的主人今日做“春日...
    開封第一講書人閱讀 31,657評論 0 22
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽泌辫。三九已至,卻和暖如春九默,著一層夾襖步出監(jiān)牢的瞬間震放,已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 32,811評論 1 268
  • 我被黑心中介騙來泰國打工驼修, 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留殿遂,地道東北人。 一個月前我還...
    沈念sama閱讀 47,685評論 2 368
  • 正文 我出身青樓乙各,卻偏偏與公主長得像墨礁,于是被迫代替她去往敵國和親。 傳聞我的和親對象是個殘疾皇子耳峦,可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 44,573評論 2 353

推薦閱讀更多精彩內(nèi)容