學(xué)校網(wǎng)課筆記~
今天的網(wǎng)課主要介紹怎么把原始的fastq文件根據(jù)fastqc結(jié)果進行去接頭、以及過濾掉低質(zhì)量的reads。
主講人介紹目前主要使用的去接頭/過濾的軟件有三種：
- TrimGalore (+ CutAdapt) https://github.com/FelixKrueger/TrimGalore/blob/master/Docs/Trim_Galore_User_Guide.md
- Trimmomatic
http://www.usadellab.org/cms/?page=trimmomatic - Fastx-toolkit
http://hannonlab.cshl.edu/fastx_toolkit/
主講人推薦的是第一個軟件。TrimGalore軟件里各種參數(shù)詳解在這里：here
先看一下原始fastq文件的QC結(jié)果：
可以看到這里有adapter是需要去掉的，有些read的質(zhì)量也不太好。還是用腳本來寫命令行（服務(wù)器運行）：
#!/bin/bash
#SBATCH --job-name=trimming_example # Job name
#SBATCH --mail-type=END,FAIL # Mail events (NONE, BEGIN, END, FAIL, ALL)
#SBATCH --mail-user=***************** # Where to send mail
#SBATCH --ntasks=1 # Run on a single CPU
#SBATCH --mem=4gb # Job memory request
#SBATCH --time=02:00:00 # Time limit hrs:min:sec
#SBATCH --output=/gpfs/home/ID/log_files/trimming_%j.log # Standard output and error log
#SBATCH -p cpu_short
module load trimgalore/0.5.0
module load python/cpu/2.7.15-ES
module load fastqc/0.11.7
# Trimming
trim_galore --paired --q 30 --gzip --fastqc /gpfs/home/ID/download/SRR1523671_1.fastq.gz /gpfs/home/ID/download/SRR1523671_2.fastq.gz -o /gpfs/home/ID/trimmed_files/
#QC again
fastqc -o QC/ /gpfs/home/ID/trimmed_files/*
上面需要注意的參數(shù)可能就是--paired和--q 30這兩個參數(shù)，因為我的文件是雙端測序的結(jié)果塘雳，所以要注明是--paired，并且要跟著兩個fastq文件普筹。還有就是你要過濾的read質(zhì)量败明，這里有個表：
一般過濾的標準是q 30 或者是q 20，一般不用40,50,60太防，如果標準設(shè)置過高妻顶，你有可能就把一些重要的信息給過濾掉了
運行后會生成一個txt文件（每一個fastq生成一個文件），是你過濾fastq后的一個summary：
SUMMARISING RUN PARAMETERS
==========================
Input filename: /gpfs/home/ID/download/SRR1523671_1.fastq.gz
Trimming mode: paired-end
Trim Galore version: 0.5.0
Cutadapt version: 1.8.2
Quality Phred score cutoff: 30
Quality encoding type selected: ASCII+33
Adapter sequence: 'AGATCGGAAGAGC' (Illumina TruSeq, Sanger iPCR; auto-detected) #這里會注明你用的是哪一種adaptor
Maximum trimming error rate: 0.1 (default)
Minimum required adapter overlap (stringency): 1 bp
Minimum required sequence length for both reads before a sequence pair gets removed: 20 bp
Running FastQC on the data once trimming has completed
Output file will be GZIP compressed
This is cutadapt 1.8.2 with Python 2.7.15 #這個軟件是要在python環(huán)境下運行的
Command line parameters: -f fastq -e 0.1 -q 30 -O 1 -a AGATCGGAAGAGC /gpfs/home/ID/download/SRR1523671_1.fas
tq.gz
Trimming 1 adapter with at most 10.0% errors in single-end mode ...
Finished in 316.18 s (15 us/read; 4.00 M reads/minute).
=== Summary ===
Total reads processed: 21,069,824 #所有的reads
Reads with adapters: 5,833,844 (27.7%) #有adaptor的reads
Reads written (passing filters): 21,069,824 (100.0%)
Total basepairs processed: 2,128,052,224 bp
Quality-trimmed: 387,450,070 bp (18.2%)
Total written (filtered): 1,719,160,219 bp (80.8%)
=== Adapter 1 ===
Sequence: AGATCGGAAGAGC; Type: regular 3'; Length: 13; Trimmed: 5833844 times.
No. of allowed errors:
0-9 bp: 0; 10-13 bp: 1
Bases preceding removed adapters:
A: 22.2%
C: 34.9%
G: 22.0%
T: 20.3%
none/other: 0.6%
Overview of removed sequences #這里是一長串的記錄杏头，有哪些reads被去掉了
length count expect max.err error counts
1 3396604 5267456.0 0 3396604
2 1146056 1316864.0 0 1146056
3 310513 329216.0 0 310513
4 109696 82304.0 0 109696
......
RUN STATISTICS FOR INPUT FILE: /gpfs/home/ID/download/SRR1523671_1.fastq.gz
=============================================
21069824 sequences processed in total
過濾后盈包，再做一次QC，做個對比：
