3 wes测序质量的控制
原视频6:测序质量的控制首先建立文件夹$ cd ~/project/wes/$ mkdir {raw,clean,align,mutation,qc}这部分包括fastqc和multiqc两个软件查看测序质量,以及使用trim_galore软件进行过滤低质量reads和去除接头。1 QC1.1 fastqc没有原视频中文件,我用了下载的三个文件做例子。乳腺癌的组织样本。所以原视频中命令我也用不上,但是还是列出来find /public/project/wes/raw_data -name *.gz|grep -v '\._'|xargs fastqc -t 10 -o ./想知道为什么要-v ._,去看原视频中的文件命名。我的文件没那么复杂,可以下面这样$ find *.gz|xargs fastqc -t 20-t 20 一次运行20个文件。如果你有很多很多文件,参考我这篇批量对多个测序文件进行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结果如下所示:
multiqc结果
Per Sequence Quality Scores
接头以上结果发现,质量可以但是需要去接头2 trim-galore 过滤低质量reads和去接头Trim Galore! is a wrapper script to automate quality and adapter trimming as well as quality control#安装conda install -c bioconda trim-galorels /path/to/your/directory/*_1.fastq.gz >1ls /path/to/your/directory/*_2.fastq.gz >2paste 1 2 > config也可以用ls|grep >打开qc.sh,写入以下内容source activate wesbin_trim_galore=trim_galoredir='/mnt/f/kelly/bioTree/server/wesproject/raw_data'cat config |while read iddoarr=(${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运行上面脚本bash qc.sh就可以了。可以看到后台正在运行kelly@DESKTOP-MRA1M1F:~$ ps -efUID PID PPID C STIME TTY TIME CMDroot 1 0 0 17:06 ? 00:00:00 /initkelly 2 1 0 17:06 tty1 00:00:01 -bashkelly 527 1 5 22:46 tty1 00:00:27 perl /home/kell/miniconda3/bin/trim_galore -q 25 --phred33 --length 36kelly 528 1 5 22:46 tty1 00:00:30 perl /home/kell/miniconda3/bin/trim_galore -q 25 --phred33 --length 36kelly 529 1 5 22:46 tty1 00:00:29 perl /home/kell/miniconda3/bin/trim_galore -q 25 --phred33 --length 36kelly 530 1 5 22:46 tty1 00:00:26 perl /home/kell/miniconda3/bin/trim_galore -q 25 --phred33 --length 36运行时间大概4h,生成以下几个文件├── [ 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.gzTrimming mode: paired-endTrim Galore version: 0.6.2Cutadapt version: 1.18Number of cores used for trimming: 1Quality Phred score cutoff: 25Quality encoding type selected: ASCII+33Adapter sequence: 'AGATCGGAAGAGC' (Illumina TruSeq, Sanger iPCR; auto-detected)Maximum trimming error rate: 0.1 (default)Minimum required adapter overlap (stringency): 3 bpMinimum required sequence length for both reads before a sequence pair gets removed: 36 bpOutput file will be GZIP compressedThis is cutadapt 1.18 with Python 2.7.15Command line parameters: -j 1 -e 0.1 -q 25 -O 3 -a AGATCGGAAGAGC SRR7696207_1.fastq.gzProcessing 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,979Reads with adapters: 9,452,510 (33.9%)Reads written (passing filters): 27,850,979 (100.0%)Total basepairs processed: 4,177,646,850 bpQuality-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: 1Bases preceding removed adapters:A: 21.7%C: 25.2%G: 28.1%T: 24.9%none/other: 0.0%Overview of removed sequenceslength count expect max.err error counts3 619264 435171.5 0 6192644 291812 108792.9 0 2918125 236788 27198.2 0 2367886 226036 6799.6 0 2260367 217969 1699.9 0 2179698 214073 425.0 0 2140739 230269 106.2 0 229591 67810 216909 26.6 1 209213 769611 223835 6.6 1 214662 917312 219837 1.7 1 210062 977513 219106 0.4 1 209421 968514 223045 0.4 1 212324 1072115 218505 0.4 1 208196 1030916 224812 0.4 1 213584 1122817 228422 0.4 1 216425 1199718 214056 0.4 1 204368 968819 216385 0.4 1 206368 1001720 207262 0.4 1 198505 875721 220284 0.4 1 209515 1076922 207937 0.4 1 198723 921423 203136 0.4 1 194604 853224 208015 0.4 1 198522 949325 200567 0.4 1 191904 866326 201338 0.4 1 192675 8663......