原視頻6:測(cè)序質(zhì)量的控制
首先建立文件夾
$ cd ~/project/wes/
$ mkdir {raw,clean,align,mutation,qc}
這部分包括fastqc和multiqc兩個(gè)軟件查看測(cè)序質(zhì)量咽安，以及使用trim_galore軟件進(jìn)行過濾低質(zhì)量reads和去除接頭。
1 QC
1.1 fastqc
沒有原視頻中文件蓬推，我用了下載的三個(gè)文件做例子妆棒。乳腺癌的組織樣本。所以原視頻中命令我也用不上沸伏，但是還是列出來
find /public/project/wes/raw_data -name *.gz|grep -v '\._'|xargs fastqc -t 10 -o ./
想知道為什么要-v ._糕珊，去看原視頻中的文件命名。
我的文件沒那么復(fù)雜毅糟，可以下面這樣
$ find *.gz|xargs fastqc -t 20
-t 20 一次運(yùn)行20個(gè)文件红选。
如果你有很多很多文件，參考我這篇批量對(duì)多個(gè)測(cè)序文件進(jìn)行fastqc.
1.2 multiqc
$ multiqc -n wes ./
......
[INFO ] multiqc : This is MultiQC v1.0.dev0
[INFO ] multiqc : Template : default
[INFO ] multiqc : Searching './'
[INFO ] fastqc : Found 8 reports
[INFO ] multiqc : Report : wes.html
[INFO ] multiqc : Data : wes_data
[INFO ] multiqc : MultiQC complete
結(jié)果如下所示：
multiqc結(jié)果
Per Sequence Quality Scores
接頭
以上結(jié)果發(fā)現(xiàn)姆另，質(zhì)量可以但是需要去接頭
2 trim-galore 過濾低質(zhì)量reads和去接頭
Trim Galore! is a wrapper script to automate quality and adapter trimming as well as quality control
#安裝
conda install -c bioconda trim-galore
ls /path/to/your/directory/*_1.fastq.gz >1
ls /path/to/your/directory/*_2.fastq.gz >2
paste 1 2 > config
也可以用
ls|grep >
打開qc.sh喇肋，寫入以下內(nèi)容
source activate wes
bin_trim_galore=trim_galore
dir='/mnt/f/kelly/bioTree/server/wesproject/raw_data'
cat config |while read id
do
arr=(${id})
fq1=${arr[0]}
fq2=${arr[1]}
nohup $bin_trim_galore -q 25 --phred33 --length 36 -e 0.1 --stringency 3 --paired -o $dir $fq1 $fq2 &
done
運(yùn)行上面腳本bash qc.sh就可以了迹辐。
可以看到后臺(tái)正在運(yùn)行
kelly@DESKTOP-MRA1M1F:~$ ps -ef
UID PID PPID C STIME TTY TIME CMD
root 1 0 0 17:06 ? 00:00:00 /init
kelly 2 1 0 17:06 tty1 00:00:01 -bash
kelly 527 1 5 22:46 tty1 00:00:27 perl /home/kelly/miniconda3/bin/trim_galore -q 25 --phred33 --length 36
kelly 528 1 5 22:46 tty1 00:00:30 perl /home/kelly/miniconda3/bin/trim_galore -q 25 --phred33 --length 36
kelly 529 1 5 22:46 tty1 00:00:29 perl /home/kelly/miniconda3/bin/trim_galore -q 25 --phred33 --length 36
kelly 530 1 5 22:46 tty1 00:00:26 perl /home/kelly/miniconda3/bin/trim_galore -q 25 --phred33 --length 36
運(yùn)行時(shí)間大概4h明吩，生成以下幾個(gè)文件
├── [ 88] 1
├── [ 88] 2
├── [ 176] config
├── [ 50K] nohup.out
├── [2.2G] SRR7696207_1.fastq.gz
├── [4.8K] SRR7696207_1.fastq.gz_trimming_report.txt
├── [2.1G] SRR7696207_1_val_1.fq.gz
├── [2.4G] SRR7696207_2.fastq.gz
├── [5.0K] SRR7696207_2.fastq.gz_trimming_report.txt
├── [2.3G] SRR7696207_2_val_2.fq.gz
├── [1.9G] SRR8517853_1.fastq.gz
├── [4.8K] SRR8517853_1.fastq.gz_trimming_report.txt
├── [1.8G] SRR8517853_1_val_1.fq.gz
├── [2.1G] SRR8517853_2.fastq.gz
├── [5.0K] SRR8517853_2.fastq.gz_trimming_report.txt
├── [2.0G] SRR8517853_2_val_2.fq.gz
├── [2.3G] SRR8517854_1.fastq.gz
├── [4.8K] SRR8517854_1.fastq.gz_trimming_report.txt
├── [2.2G] SRR8517854_1_val_1.fq.gz
├── [2.6G] SRR8517854_2.fastq.gz
├── [5.1K] SRR8517854_2.fastq.gz_trimming_report.txt
├── [2.4G] SRR8517854_2_val_2.fq.gz
├── [4.1G] SRR8517856_1.fastq.gz
├── [4.9K] SRR8517856_1.fastq.gz_trimming_report.txt
├── [4.0G] SRR8517856_1_val_1.fq.gz
├── [4.5G] SRR8517856_2.fastq.gz
├── [5.1K] SRR8517856_2.fastq.gz_trimming_report.txt
├── [4.2G] SRR8517856_2_val_2.fq.gz
└── [ 305] trim
可見贺喝，各小了大概0.1G。
其中氮采，txt中的信息如下
SUMMARISING RUN PARAMETERS
==========================
Input filename: SRR7696207_1.fastq.gz
Trimming mode: paired-end
Trim Galore version: 0.6.2
Cutadapt version: 1.18
Number of cores used for trimming: 1
Quality Phred score cutoff: 25
Quality encoding type selected: ASCII+33
Adapter sequence: 'AGATCGGAAGAGC' (Illumina TruSeq, Sanger iPCR; auto-detected)
Maximum trimming error rate: 0.1 (default)
Minimum required adapter overlap (stringency): 3 bp
Minimum required sequence length for both reads before a sequence pair gets removed: 36 bp
Output file will be GZIP compressed
This is cutadapt 1.18 with Python 2.7.15
Command line parameters: -j 1 -e 0.1 -q 25 -O 3 -a AGATCGGAAGAGC SRR7696207_1.fastq.gz
Processing reads on 1 core in single-end mode ...
Finished in 1011.58 s (36 us/read; 1.65 M reads/minute).
=== Summary ===
Total reads processed: 27,850,979
Reads with adapters: 9,452,510 (33.9%)
Reads written (passing filters): 27,850,979 (100.0%)
Total basepairs processed: 4,177,646,850 bp
Quality-trimmed: 21,999,885 bp (0.5%)
Total written (filtered): 3,930,084,164 bp (94.1%)
=== Adapter 1 ===
Sequence: AGATCGGAAGAGC; Type: regular 3'; Length: 13; Trimmed: 9452510 times.
No. of allowed errors:
0-9 bp: 0; 10-13 bp: 1
Bases preceding removed adapters:
A: 21.7%
C: 25.2%
G: 28.1%
T: 24.9%
none/other: 0.0%
Overview of removed sequences
length count expect max.err error counts
3 619264 435171.5 0 619264
4 291812 108792.9 0 291812
5 236788 27198.2 0 236788
6 226036 6799.6 0 226036
7 217969 1699.9 0 217969
8 214073 425.0 0 214073
9 230269 106.2 0 229591 678
10 216909 26.6 1 209213 7696
11 223835 6.6 1 214662 9173
12 219837 1.7 1 210062 9775
13 219106 0.4 1 209421 9685
14 223045 0.4 1 212324 10721
15 218505 0.4 1 208196 10309
16 224812 0.4 1 213584 11228
17 228422 0.4 1 216425 11997
18 214056 0.4 1 204368 9688
19 216385 0.4 1 206368 10017
20 207262 0.4 1 198505 8757
21 220284 0.4 1 209515 10769
22 207937 0.4 1 198723 9214
23 203136 0.4 1 194604 8532
24 208015 0.4 1 198522 9493
25 200567 0.4 1 191904 8663
26 201338 0.4 1 192675 8663
......
后面如果想快速運(yùn)行流程主到，可以把測(cè)序數(shù)據(jù)取前N行躯概，那么請(qǐng)看
3_0_4 要理解并會(huì)用的幾個(gè)腳本
如果就想運(yùn)行所有數(shù)據(jù)，請(qǐng)到
4 比對(duì)到參考基因組輸出bam文件
