DAY5學(xué)習(xí)內(nèi)容：過濾和比對
過濾
通過DAY4_fastqc結(jié)果裳食，我們開始對數(shù)據(jù)進(jìn)行過濾茧球，主要是去除接頭和低質(zhì)量堿基蝉绷，一般使用兩種軟件可以完成這個(gè)過程：trim_galore和trimmoatic拜姿，而這兩種軟件遵循的原則是直接從不合格的位置向后全部切掉驻粟，而不是只挖掉這塊然评，為了避免人為改變測序序列尘惧。
兩款軟件的學(xué)習(xí)資料：
trim_galore：https://github.com/FelixKrueger/TrimGalore/blob/master/Docs/Trim_Galore_User_Guide.md
trimmomatic：https://github.com/mossmatters/KewHybSeqWorkshop/blob/master/FastQC_Trimmomatic.md
trim_galore
trim_galore在去除接頭序列時(shí)依賴的是Cutadapt炬转，一般接頭出現(xiàn)在3'末端；
選項(xiàng)參數(shù)：
-q：設(shè)置堿基質(zhì)量閾值
--phredxx：質(zhì)量體系述吸；phred33和phred64忿族，不過現(xiàn)在絕大部分 Illumina 平臺的產(chǎn)出數(shù)據(jù)也都轉(zhuǎn)為使用 phred33 格式了。
如何區(qū)分phred33和phred64：fa第四行的堿基質(zhì)量值如果大部分都是大寫字母那就是phred64刚梭，如果存在特殊符號類似#肠阱、%票唆、@就是phred33朴读；或者看fastqc報(bào)告中basic statistics結(jié)果：Illumina 1.8版本以上的就是phred33
--length：表示測序片段長度閾值，例如--length 50走趋，表示如果測序長度為150bp衅金，如果左右切掉接頭和低質(zhì)量堿基后，長度低于了50bp簿煌，就直接把這整條去掉
--stringency：設(shè)置的數(shù)字表示：認(rèn)為有幾個(gè)堿基和接頭序列是有重疊的氮唯，默認(rèn)值為1，意思是從3'的對應(yīng)位置檢測姨伟，出現(xiàn)一個(gè)堿基與常用接頭有重疊就把這個(gè)堿基以后的序列都去掉惩琉。
常用的接頭前12-13bp序列
Illumina: AGATCGGAAGAGC
Small RNA: TGGAATTCTCGG
Nextera: CTGTCTCTTATA
--paired：雙末端測序
-o：輸出路徑
以上參數(shù)及其他參數(shù)執(zhí)行trim_galore -help均可查看
trim_galore -help
通過幫助文件了解軟件使用后，我們開始用起它夺荒，兩種情況：
-第一種情況：只對一組fq文件進(jìn)行過濾
$dir=/YOUR_PATH/ #定義輸出文件路徑
$fq1=/YOUR_PATH/xxx.fq1 #將xxx.fq1文件路徑及本身賦值給fq1
$fq2=/YOUR_PATH/xxx.fq2 #將xxx.fq2文件路徑及本身賦值給fq2
trim_galore -q 25 --phred33 --length 50 -e 0.1 --stringency 3 --paired -o $dir $fq1 $fq2 #-e 允許的最大錯(cuò)誤率（錯(cuò)誤數(shù)除以匹配區(qū)域的長度）（默認(rèn)值：0.1）
-第二種情況：批量對多組fq文件進(jìn)行過濾
首先要將所有的fq文件放在同一個(gè)文件中瞒渠，基于前期學(xué)習(xí)我們這共有8組fq文件，如下：
如何快速的完成這個(gè)過程呢技扼？
#進(jìn)入fq的文件目錄
cd raw
#找到1端測序文件伍玖，并保存到另一個(gè)文件中，例如這里叫做fastq_1
ls *_1* >fastq_1
#找到2端測序文件剿吻，并保存到另一個(gè)文件中窍箍，例如這里叫做fastq_2
ls *_2* >fastq_2
#接下來我們將這兩個(gè)新生成的文件合并起來，使用paste命令
paste fastq_1 fastq_2
#屏幕顯示所有fq文件
為了保證每個(gè)文件在任何時(shí)候和目錄下都可以重復(fù)使用，我們需要添加其絕對路徑椰棘，因此我們要對上面的工作再加上其路徑纺棺，并將其合并結(jié)果放到一個(gè)新的文件中，命令如下：
ls `pwd`/*_1* >new_fastq_1
ls `pwd`/*_2* >new_fastq_2
paste new_fastq_1 new_fastq_2 >conf_fastq
#conf意思是配置文件的意思邪狞，然后#cat conf_fastq五辽，顯示加上了路徑
對上面的conf_fastq進(jìn)行trim_galore操作，同樣使用循環(huán)命令外恕，如下：
#將腳本寫入script目錄中杆逗，
cd script
#新建一個(gè)名字為filter.sh腳本，并添加內(nèi)容
cat >filter.sh #或
vi filter.sh
#腳本內(nèi)容
filter_data=~/raw
cat $filter_data/conf_fastq | while read i
do
fqs=($i)
fq1=${fqs[0]}
fq2=${fqs[1]}
echo $i
# nohup trim_galore -q 25 --phred33 --length 50 -e 0.1 --stringency 3 --paired -o /hwfssz1/ST_AGRIC/F13ZHQYJSY1409_LP/USER/liujia/TCM/TD/clean $fq1 $fq2 &
done
還是先將主要執(zhí)行的命令注釋掉鳞疲，看看循環(huán)可行不罪郊，如下：
沒問題！釋放執(zhí)行腳本尚洽！
filter_data=~/raw
cat $filter_data/conf_fastq | while read i
do
fqs=($i)
fq1=${fqs[0]}
fq2=${fqs[1]}
# echo $i
nohup trim_galore -q 25 --phred33 --length 50 -e 0.1 --stringency 3 --paired -o /hwfssz1/ST_AGRIC/F13ZHQYJSY1409_LP/USER/liujia/TCM/TD/clean $fq1 $fq2 &
done
在上述腳本中悔橄，一個(gè)新學(xué)的只是就是創(chuàng)建數(shù)組，進(jìn)行數(shù)據(jù)提取
fqs=($i)
fq1=${fqs[0]}
fq2=${fqs[1]}
#使用（）創(chuàng)建數(shù)組腺毫，將每一個(gè)$i（一對PE測序文件）放入其中
#然后使用${[]}進(jìn)行訪問癣疟，[]表示訪問第幾列，從0開始計(jì)數(shù)
#每當(dāng)讀入一對fq時(shí)潮酒，${fqs[0]}表示將第一個(gè)測序文件也就是_1文件賦值給fq1睛挚，同樣${fqs[1]}表示將第二個(gè)測序文件也就是_2文件賦值給fq2，此過程必須要帶絕對路徑急黎。
腳本編輯完扎狱，我們就可以運(yùn)行腳本了，bash filter.sh 或 sh filter.sh
sh和bash區(qū)別：
sh勃教、bash都是linux中的解釋器淤击，但有的sh是沒有數(shù)組array解釋器的，因此當(dāng)你使用sh去運(yùn)行整個(gè)命令時(shí)故源，有可能會產(chǎn)生報(bào)錯(cuò)：
Syntax error: "(" unexpected (expecting "done")
但這并不是說代碼寫錯(cuò)了污抬，而是選擇錯(cuò)了解釋器。
trimmomatic：
該軟件有兩種過濾模式绳军，分別對應(yīng)SE和PE測序數(shù)據(jù)印机，同時(shí)支持gzip和bzip2壓縮文件。Trimmomatic 過濾數(shù)據(jù)的步驟與命令行中過濾參數(shù)的順序有關(guān)
過濾步驟：
- ILLUMINACLIP删铃；
- SLIDINGWINDOW耳贬；
- MAXINFO；
- LEADING猎唁；
- TRAILING咒劲；
- CROP顷蟆；
- HEADCROP；
- MINLEN腐魂；
- AVGQUAL帐偎；
- TOPHRED33；
- TOPHRED64
參數(shù)選擇：
對于雙末端測序模式：
PE模式
輸入輸出文件：
PE 模式的兩個(gè)輸入文件：sample_R1.fastq sample_R2.fastq以及四個(gè)輸出文件：sample_paired_R1.clean.fastq sample_unpaired_R1.clean.fastq sample_paired_R1.clean.fastq sample_unpaired_R1.clean.fastq
SLIDINGWINDOW：滑窗剪切蛔屹，統(tǒng)計(jì)滑窗口中所有堿基的平均質(zhì)量值削樊，如果低于設(shè)定的閾值，則切掉窗口兔毒。
SLIDINGWINDOW:<windowSize>:<requiredQuality>
widowSize:設(shè)置窗口大小
requiredQuality：設(shè)置窗口堿基平均質(zhì)量閾值
SLIDINGWINDOW:5:20 設(shè)置 5bp 窗口漫贞，堿基平均質(zhì)量值閾值 20
LEADING：從 reads 的起始端開始切除質(zhì)量值低于設(shè)定的閾值的堿基，直到有一個(gè)堿基其質(zhì)量值達(dá)到閾值育叁。
LEADING:<quality>
quality:設(shè)定堿基質(zhì)量值閾值迅脐，低于這個(gè)閾值將被切除。
LEADING:5
TRAILING：從 reads 的末端開始切除質(zhì)量值低于設(shè)定閾值的堿基豪嗽，直到有一個(gè)堿基質(zhì)量值達(dá)到閾值谴蔑。Illumina 平臺有些低質(zhì)量的堿基質(zhì)量值被標(biāo)記為 2 ，所以設(shè)置為 3 可以過濾掉這部分低質(zhì)量堿基龟梦。TRAILING:<quality>
quality:設(shè)定堿基質(zhì)量值閾值隐锭，低于這個(gè)閾值將被切除。
TRAILING:5
MINLEN：設(shè)定一個(gè)最短 read 長度计贰，當(dāng) reads 經(jīng)過前面的過濾之后钦睡，如果保留下來的長度低于這個(gè)閾值，整條 reads 被丟棄蹦玫。被丟棄的 reads 數(shù)會被統(tǒng)計(jì)在 Trimmomatic 日志的 dropped reads 中赎婚。
MINLEN:<length>
length：可被保留的最短 read 長度
MINLEN:20
-trimlog：輸出Trimmomatic 日志
TOPHRED33/64：
此選項(xiàng)可以將過濾之后的 Fastq 文件中質(zhì)量值這一行轉(zhuǎn)為 phred-33 格式 或 phred-64 格式或。
更多關(guān)于Trimmomatic的學(xué)習(xí)資料請看考NGS 數(shù)據(jù)過濾之 Trimmomatic 詳細(xì)說明或其--help
了解Trimmomatic軟件使用后樱溉，我們開始用起它，腳本如下纬凤，參數(shù)如上：
mkdir $rna/clean/trimmomatic && cd "$_"
cat >filter-trim.sh
rna=/YOUR_PATH/rnaseq
cat $rna/raw/conf_fastq | while read i
do
fqs=($i)
fq1=${fqs[0]}
fq2=${fqs[1]}
tri1=`basename $fq1`
tri2=`basename $fq2`
nohup trimmomatic PE -phred33 \
-trimlog trim.logfile\
$fq1 $fq2 \
clean.$tri1 unpaired.$tri1 \
clean.$tri2 unpaired.$tri2 \
SLIDINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:20 &
done
bash filter-trim.sh
任務(wù)已經(jīng)提交福贞，那我們怎么看我們這個(gè)任務(wù)跑的是否成功呢
利用htop工具（進(jìn)程管理器）功能和top一樣，只是它是彩色的停士，信息詳細(xì)
安裝學(xué)習(xí)鏈接：
conda安裝看這里：https://anaconda.org/conda-forge/htop
官網(wǎng)安裝看這里：https://hisham.hm/htop/
安裝后挖帘，直接輸入htop回車就可以查看當(dāng)前進(jìn)程狀態(tài)：CPU、內(nèi)存恋技、正在跑的進(jìn)程
如果你只想看你自己的進(jìn)程狀態(tài)：htop -u your_ID
過濾好的數(shù)據(jù)我們還要進(jìn)行質(zhì)控拇舀，看看接頭和低質(zhì)量的堿基是否被處理掉，主要看接頭蜻底！
比對
基于之前工作得到的cleandata骄崩，開始進(jìn)行比對過程
比對所需的文件：
-過濾并質(zhì)控合格后的數(shù)據(jù)
-參考基因組和注釋文件
如何生成小數(shù)據(jù)用于測試流程
生成小數(shù)據(jù)，說白了就是從你的數(shù)據(jù)里提取出一部分，用于測試要拂，這里我們使用seqtk軟件抠璃，https://github.com/lh3/seqtk，用conda安裝即可脱惰，
首先我們先整理一下之前生成的cleandata搏嗡，方法如上：
# 先配置cleandata數(shù)據(jù)集文件
cd clean
ls `pwd`/*_1*gz >fastqc_1.clean
ls `pwd`/*_2*gz >fastqc_2.clean
paste fastqc_1.clean fastqc_2.clean >fastq_clean.conf
cat fastq_clean.conf
在項(xiàng)目總的目錄下新建一個(gè)test目錄
#創(chuàng)建test目錄并進(jìn)入
mkdir test && cd $_
#vi編譯名字為sample_fq.sh腳本
vi sample_fq.sh
#腳本內(nèi)容
clean_data=~/clean
cat $clean_data/fastq_clean.conf | while read i;
do
fqs=($i)
fq1=${fqs[0]}
fq2=${fqs[1]}
file=`basename $fq1`
#echo $file
surname=${file%%_*}
#echo $fq1 $fq2 $surname
# 隨機(jī)選了10000條reads
seqtk sample $fq1 10000 >test.${surname}_1.fastq
seqtk sample $fq2 10000 >test.${surname}_2.fastq
done
bash sample_fq.sh
以上腳本可以通過echo $file，可得知file做了什么拉一，basename就是數(shù)據(jù)這個(gè)路徑下最底層的名字采盒，輸出如下：
echo $surname又做了什么？輸出如下
$fq1：~/clean/SRR1039508_1_val_1. fq.gz
$fq2：~/clean/SRR1039508_2_val_2. fq.gz
$surname：SRR1039508
file%%_*：%%就是將file中_符號后面的全部刪除蔚润，也就是把_2_val_2. fq.gz刪除纽甘，只留下SRR1039508
生成小數(shù)據(jù)結(jié)果如下：
比對第一課：hisat2比對
比對又叫mapping，mapping的目的就是找到reads在基因組上的位置抽碌。
hisat2官網(wǎng)：http://ccb.jhu.edu/software/hisat2/index.shtml
選項(xiàng)參數(shù)：http://ccb.jhu.edu/software/hisat2/manual.shtml
文章推薦：NC的文章：Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis
更多背景知識請自行查找悍赢，開始對上一步生成的小數(shù)據(jù)進(jìn)行mapping
mapping第一步：構(gòu)建基因組索引
構(gòu)建索引的目的：快速在基因組上進(jìn)行定位，因?yàn)閞eads一般都是幾千萬條货徙，基因組也有30億堿基左权，肯定不能一條reads在基因組遍歷一遍，再進(jìn)行下一條reads痴颊。
mkdir ~/align/hisat2 && cd "$_"
vi index.sh #創(chuàng)建index.sh腳本
#腳本內(nèi)容
rna=/YOUR_PATH/rnaseq
genome=$rna/ref/hg19.fa #DAY4下載的參考基因組
hisat2-build -p 10 $genome hg19 #建立索引
mapping第二步：開始比對
建立好索引后赏迟，我們開始進(jìn)行比對（其實(shí)第二步在mapping之后還包含sam轉(zhuǎn)bam和排序），腳本如下蠢棱，改腳本在~/align/hisat2路徑下執(zhí)行：
# 小數(shù)據(jù)集合的路徑
ls ~/test/*_1* >hisat2_1.test
ls ~/test/*_2* >hisat2_2.test
paste 1.test 2.test >test_hisat2.conf
vi hisat2_aln.sh
cat ~/test/test_hisat2.conf
| while read i
do
fqs=($i)
fq1=${fqs[0]}
fq2=${fqs[1]}
sample=`basename $fq1 | cut -d '_' -f1`
hisat2 -p 10 -x hg19 -1 $fq1 -2 $fq2 | samtools sort -T PREFIX.nnnn.bam -O bam \
-@ 10 -o ${sample}_hisat.sorted.bam
samtools index -b ${sample}_hisat.sorted.bam
done
參數(shù)信息：
-p：進(jìn)程數(shù)锌杀，越打越快
-x：參考基因組索引文件的前綴
-1：雙端測序結(jié)果的第一個(gè)文件。若有多組數(shù)據(jù)泻仙，使用逗號將文件分隔糕再。Reads的長度可以不一致。
-2：雙端測序結(jié)果的第二個(gè)文件玉转。若有多組數(shù)據(jù)突想，使用逗號將文件分隔，并且文件順序要和-1參數(shù)對應(yīng)究抓。Reads的長度可以不一致猾担。
samtools sort：使用samtools工具對生成的sam文件進(jìn)行排序，并轉(zhuǎn)換成二進(jìn)制bam文件
-T：添加一個(gè)臨時(shí)文件是必須的刺下，格式：PREFIX.nnnn.bam
-O：設(shè)置輸出格式包括：bam绑嘹，sam或cram
-@：除主線程外還要使用的BAM壓縮線程數(shù) 10
samtools index：對bam進(jìn)行排好隊(duì)后，需要給大家一個(gè)編號橘茉，就是index
-b：建立索引后將產(chǎn)生后綴為.bai的文件工腋。
以上腳本工作流程如下：hisat2將reads比對到了參考基因組姨丈，理論上會生成sam文件，但是這個(gè)文件太大夷蚊，因此通過samtools sort -O bam直接轉(zhuǎn)換成二進(jìn)制的bam文件并排序构挤，與sam內(nèi)容一樣而且節(jié)省空間。最后為了下一步導(dǎo)入進(jìn)行可視化并且為了定量的需要惕鼓，我們又對bam構(gòu)建了索引（samtools index）筋现。
對于samtools安裝，conda install samtools=1.8箱歧，一定要設(shè)置1.8安裝版本矾飞，不然會報(bào)錯(cuò)
關(guān)于更多samtools請學(xué)習(xí)：處理SAM、BAM你需要Samtools呀邢；英文資料：http://www.htslib.org/doc/
比對后生成文件洒沦，mapping率在nohup.out下看（less nohup.out）如下：
