這里是佳奧利凑！繼續(xù)轉(zhuǎn)錄組分析的學(xué)習(xí)。
一般我們都是從.fq文件開始的，但是如果是自己尋找數(shù)據(jù)的話哀澈，要多一步驟牌借。
1 獲取fastq數(shù)據(jù)
1.1 sratoolkit
我們看文章的時(shí)候，找到對(duì)應(yīng)的SRR號(hào)割按。
進(jìn)入環(huán)境走哺，下載sratoolkit軟件，后續(xù)就可以直接在Linux內(nèi)操作了哲虾。
最好的方式還是從官網(wǎng)下載（Github）丙躏，把壓縮包傳到Linux，解壓（conda下載的無法使用束凑，版本過舊）
添加文件夾到環(huán)境變量即可
$ export PATH="$PATH:/home/kaoku/biosoft/sratoolkit/sratoolkit.3.0.0-ubuntu64/bin"
配置軟件：這個(gè)界面是可以鼠標(biāo)點(diǎn)擊的晒旅，設(shè)置兩個(gè)路徑到root/ncbi即可
$ vdb-config --interactive
安裝成功后，使用prefetch便開始下載
prefetch SRR1039508
prefetch批量下載：SRR存在一個(gè)id.txt內(nèi)
cat id | while read id; do prefetch $id; done
掛在后臺(tái)下載：
cat id | while read id; do (prefetch $id &); done
使用top查看下載進(jìn)程汪诉。
想關(guān)掉全部進(jìn)程：
ps -ef | grep prefe | awk '{print $2}'| while read id; do kill $id; done
下載完成后废恋，我們會(huì)看到一列的SRR1039508.sra文件，轉(zhuǎn)化成fastq文件
fastq-dump --gzip --split-3 -o ./ /SRR1039508.sra
便開始生成SRR1039508_1.fastq.gz文件扒寄。
1.2 直接wget
這是一個(gè)研究對(duì)象為擬南芥的文章鱼鼓，所有的fastq數(shù)據(jù)存放于此， ID為E-MTAB-5130该编。
先獲取.txt文件迄本，再提取出URL，wget下載课竣。
wget http://www.ebi.ac.uk/arrayexpress/files/E-MTAB-5130/E-MTAB-5130.sdrf.txt
head -n1 E-MTAB-5130.sdrf.txt | tr '\t' '\n' | nl | grep URI
tail -n +2 E-MTAB-5130.sdrf.txt | cut -f 33 | xargs -i wget {}
隨后下載參考基因組TAIR10 ver.24嘉赎。
nohup wget ftp://ftp.ensemblgenomes.org/pub/plants/release-28/fasta/arabidopsis_thaliana/cdna/Arabidopsis_thaliana.TAIR10.28.cdna.all.fa.gz &
nohup wget ftp://ftp.ensemblgenomes.org/pub/plants/release-28/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.28.dna.genome.fa.gz &
nohup wget ftp://ftp.ensemblgenomes.org/pub/plants/release-28/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.28.gff3.gz &
nohup wget ftp://ftp.ensemblgenomes.org/pub/plants/release-28/gtf/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.28.gtf.gz &
后續(xù)就是按照流程一步一步了。
查看下載日志：
less raw/nohup.out
grep failed raw/nohup.out
2 fastqc查看測(cè)序質(zhì)量
獲取數(shù)據(jù)后于樟，用fastqc查看數(shù)據(jù)質(zhì)量公条。
fastqc SRR957678.fastq.gz
批量處理.fastqc.gz：一共10個(gè)文件
ls *gz | xargs fastqc -t 10
用multiqu批量查看.html結(jié)果
使用conda安裝：
conda install -c bioconda multiqc
進(jìn)入結(jié)果目錄，直接運(yùn)行：
multiqc ./
即可出現(xiàn)multiqc_report.html文件
下面我用數(shù)據(jù)跑一遍流程：
$ ls
SRR957677.fastq.gz SRR957678.fastq.gz SRR957679.fastq.gz SRR957680.fastq.gz
全部跑一遍fastqc
$ ls *gz | xargs fastqc -t 10
運(yùn)行multqc
$ multiqc ./
/// MultiQC | v1.13.dev0
| multiqc | Search path : /home/kaoku/rnaseq/biotree_human
| searching | ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 14/14
| fastqc | Found 4 reports
| multiqc | Compressing plot data
| multiqc | Report : multiqc_report.html
| multiqc | Data : multiqc_data
| multiqc | MultiQC complete
$ ls 查看該multiqc_report.html
multiqc_data
multiqc_report.html
SRR957677_fastqc.html
SRR957677_fastqc.zip
SRR957677.fastq.gz
SRR957678_fastqc.html
SRR957678_fastqc.zip
SRR957678.fastq.gz
SRR957679_fastqc.html
SRR957679_fastqc.zip
SRR957679.fastq.gz
SRR957680_fastqc.html
SRR957680_fastqc.zip
SRR957680.fastq.gz
整體質(zhì)量還不錯(cuò)迂曲。
查看了數(shù)據(jù)以后靶橱，我們就需要過濾測(cè)序數(shù)據(jù)。
3 過濾測(cè)序結(jié)果（雙端測(cè)序結(jié)果）
如果Adapter Content是錯(cuò)誤的路捧，說明需要對(duì)reads進(jìn)行修剪處理关霸。
查看Adapter Content ：（fastqc軟件測(cè)出的Adapter Content 序列后數(shù)個(gè)堿基是TCTTCTGCTTG）
zcat SRR957677_fastqc.zip | grep TCTTCTGCTTG
使用TrimGalore軟件
在GitHub下載.tar.gz，解壓
添加到環(huán)境變量（要進(jìn)入rnaseq環(huán)境鬓长，依賴cutadapt）
$ export PATH="$PATH:/home/kaoku/biosoft/trim_galory/TrimGalore-0.6.7"
設(shè)置參數(shù)：（--length 30-40谒拴、--stringency 3-5）
dir='/home/kaoku/rnaseq/biotree_human/clean'（工作目錄）
fq1='/home/kaoku/rnaseq/biotree_human/SRR957677_1_fastqc.zip'（舉例）
fq2='/home/kaoku/rnaseq/biotree_human/SRR957677_2_fastqc.zip'（舉例）
trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 3 --paired -o $dir $fq1 $fq2
可以繼續(xù)調(diào)整參數(shù)（引自lncRNA組裝流程的軟件介紹之trim-galore）
--quality：設(shè)定Phred quality score閾值尝江，默認(rèn)為20涉波。分析時(shí)可改成25，稍微嚴(yán)格一些。
--phred33：：選擇-phred33或者-phred64啤覆，表示測(cè)序平臺(tái)使用的Phred quality score苍日。具體怎么選擇，看你用什么測(cè)序平臺(tái)窗声；例如：-phred33對(duì)應(yīng)(Sanger/Illumina 1.9+ encoding)相恃，-phred64對(duì)應(yīng)(Illumina 1.5 encoding)
--adapter：輸入adapter序列。也可以不輸入笨觅，Trim Galore!會(huì)自動(dòng)尋找可能性最高的平臺(tái)對(duì)應(yīng)的adapter拦耐。自動(dòng)搜選的平臺(tái)三個(gè)，也直接顯式輸入這三種平臺(tái)见剩，即--illumina杀糯、--nextera和--small_rna。
--stringency：設(shè)定可以忍受的前后adapter重疊的堿基數(shù)苍苞，默認(rèn)為1（非彻毯玻苛刻）「牵可以適度放寬骂际，因?yàn)楹笠粋€(gè)adapter幾乎不可能被測(cè)序儀讀到。
--length：設(shè)定輸出reads長(zhǎng)度閾值冈欢，小于設(shè)定值會(huì)被拋棄歉铝。
--paired：對(duì)于雙端測(cè)序結(jié)果，一對(duì)reads中凑耻，如果有一個(gè)被剔除犯戏，那么另一個(gè)會(huì)被同樣拋棄，而不管是否達(dá)到標(biāo)準(zhǔn)拳话。
--retain_unpaired：對(duì)于雙端測(cè)序結(jié)果先匪，一對(duì)reads中，如果一個(gè)read達(dá)到標(biāo)準(zhǔn)弃衍，但是對(duì)應(yīng)的另一個(gè)要被拋棄呀非，達(dá)到標(biāo)準(zhǔn)的read會(huì)被單獨(dú)保存為一個(gè)文件。
--gzip和--dont_gzip：清洗后的數(shù)據(jù)zip打包或者不打包镜盯。
--output_dir：輸入目錄岸裙。需要提前建立目錄，否則運(yùn)行會(huì)報(bào)錯(cuò)速缆。
--trim-n : 移除read一端的reads
批量處理的話：
ls *_1.fastqc.gz >1
ls *_2.fastqc.gz >2
paste 1 2 > config
dir='/home/kaoku/rnaseq/biotree_human/clean'
cat config | while read id
do
arr=($id)
fq1=${arr[0]}
fq2=${arr[1]}
echo $dir $fq1 $fq2
nohub trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 3 --paired -o $dir $fq1 $fq2 &
done
運(yùn)行完成后也會(huì)有ERR1698194_2_val_2_fastqc.html報(bào)告降允，查看清洗前后差異。
拿到干凈的測(cè)序數(shù)據(jù)后艺糜，就可以開始后面的比對(duì)分析流程了剧董。
那么我們幢尚，下篇再見！
