生信基础 质控的那些事

质控的方向

就平常的生物实验而言，质控是要贯穿到实验的各个阶段的，比如那些阳性阴性对照就可以认为是在进行质控了。而在测序数据的分析过程中，也同样需要将质控贯穿到分析的各个阶段。质控基本上可以分为三个阶段，依据你的分析对象和分析目的不同而有所变化：raw sequence data, alignment, variant calling. 目前比较常见的质控多是局限在raw data的质控这个阶段。虽然如此，另外两个阶段的质控也应该是需要的，因为后续的许多分析也是基于这两个阶段。Raw data 质控的方向通常有以下几类：quality trimming, adapter removal, contaminant filtering. 而在alignment阶段，我们可以通过uniquely mapped reads, signed noise ratio 进行质控。在variant calling阶段，还在调查中……

QC of Raw data

FastQC质量统计

1
fastqc /path/*.fastq.gz -o /path/QC_result/

质控条件

  • 去除非测序序列：adapter、primer、barcode、index
  • 去除低质量reads：低于Q20、Q30的碱基占比达30%
  • 去除含N过多reads: 占比达10%
  • 去除重复序列：RNA-seq、16s-seq一般不做去除
  • 去除3’端质量Q低于10的碱基，即碱基错误率为0.1
  • 去除较短长度reads：低于25的去除
  • 根据分析结果, 适当地重复上面某些步骤

使用Trimmomatic修剪reads

1
2
3
4
5
6
7
8
trimmomatic PE read1.fq.gz  read2.fq.gz \
    clean_paired_forward.fq.gz clean_unpaired_forward.fq.gz \
    clean_paired_reverse.fq.gz clean_unpaired_reverse.fq.gz \
    ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
    LEADING:4 TRAILING:5 SLIDINGWINDOW:4:5 MINLEN:25 &
# ILLUMINACLIP:adpater存在.fa文件中，允许2个错配；回文模式下匹配阈值;；简单模式下匹配阈值
# LEADING：切除首端质量值小于4的碱基
# TRAILING: 切除末端质量小于5的碱基

使用cutadapt 修剪 reads

1
cutadapt -a ADAPTER_FWD -A ADAPTER_REV -o out.1.fastq -p out.2.fastq reads.1.fastq reads.2.fastq

fastx_toolkit 个性化清洗和过滤

1
2
3
fastq_quality_filter -q 30 -p 80 -z -i input.fastq -o out.fastq
# input.fastq不支持gz，可以如下,使用-c保留源文件
gunzip -c input.fastq.gz | fastq_quality_filter ....etc....

multiqc

fastp

据说是一步到位

另外两个阶段的质控

在比对阶段和变异检测阶段的质控软件不多，已有的软件通常都包含了这两个阶段的质控，所以合并起来。

使用RSeQC进行质控

http://rseqc.sourceforge.net/#usage-information
该软件的input格式有SAM、BAM、Fasta、BED

1 2 3 4
bam2fq.py bam2wig.py bam_stat.py insertion_profile.py
clipping_profile.py deletion_profile.py divide_bam.py infer_experiment.py
FPKM_count.py eneBody_coverage.py geneBody_coverage2.py inner_distance.py
junction_annotation.py junction_saturation.py mismatch_profile.py overlay_bigwig.py
normalize_bigwig.py read_distribution.py read_duplication.py read_GC.pyread_hexamer.py
read_NVC.py read_quality.py RNA_fragment_size.py RPKM_count.py
RPKM_saturation.py spilt_bam.py split_paired_bam.py none 
1
2
3
python scripts/bam2fq.py -i test_PairedEnd_StrandSpecific_hg19.sam  -o bam2fq_out1

python clipping_profile.py -i Pairend_StrandSpecific_51mer_Human_hg19.bam -s "PE" -o out

参考

Comments

▲