第2篇:原始数据的质控、比对和过滤

专题目录:

这部分内容包括对原始测序数据质控,然后比对过滤,这是所有NGS数据处理的上游分析。

  • ATAC-Seq与其他方法不同的一点是需要过滤去除线粒体(如果是植物,还需要过滤叶绿体),因为线粒体DNA是裸露的,也可以被Tn5酶识别切割。
  • 另外一点需要注意的是课程中给出的是单端比对的示例代码,如果是双端测序做相应更改即可。

学习目标

  • 用FastQC进行质控检测
  • 用Trimmomatic进行质量过滤
  • 用Bowtie2比对,并理解相关参数含义

测序reads 的质控流程示意图

img

FASTQC

首先对拿到的原始测序数据(fastq或fastq.gz格式)进行质控检测,直接用fastqc软件,再加上multiqc将多个检测结果一起展示。 如:

fastqc -o out_dir raw_data/*gz
multiqc *fastqc.zip --ignore *.html

Trimmomatic

Trimmomatic 可以用于去除接头,过滤低质量数据。相同功能的软件还有很多,如trim_galorecutadapt等,个人比较喜欢trim_galore可以自动识别接头类型。

# 课程中给出的Trimmomatic 的用法(单端测序)
$ java -jar /opt/Trimmomatic-0.33/trimmomatic-0.33.jar SE \
-threads 2 \
-phred33 \
H1hesc_Input_Rep1_chr12.fastq \
../results/trimmed/H1hesc_Input_Rep1_chr12.qualtrim20.minlen36.fq \
LEADING:20 \
TRAILING:20 \
MINLEN:36

Trimmomatic参数含义:可以参考NGS 数据过滤之 Trimmomatic 详细说明 trim_galore使用示例

trim_galore -q 20 --phred33 --stringency 3 --length 20 -e 0.1 --paired fq1 fq2  --gzip -o input_data_dir
# 重新用fastqc检测进行过滤后的reads质量
fastqc -o out_dir *fq.gz
multiqc *fastqc.zip --ignore *.html

比对

Bowtie2是一个快速精确的比对工具,基于Burrows-Wheeler Transform 构建基因组的FM 索引,比对过程所耗内存少。Bowtie2支持局部、双端、缺口比对模式,对大于50bp的reads比对效果更好(小于50bp的reads用Bowtie1)。

创建Bowtie2索引

bowtie2-build <path_to_reference_genome.fa> <prefix_to_name_indexes>
# Can find indexes for the entire genome on Orchestra using following path: /groups/shared_databases/igenome/Homo_sapiens/UCSC/hg19/Sequence/Bowtie2Index/

Bowtie2 比对

  • p: 线程数
  • q: reads是fastq格式
  • x: index路径
  • U: fastq路径
  • S: 输出Sam格式文件
## 课程中给出的代码是单端比对
bowtie2 -p 2 -q \
-x ~/ngs_course/chipseq/reference_data/chr12\
-U ~/ngs_course/chipseq/results/trimmed/H1hesc_Input_Rep1_chr12.qualtrim20.minlen36.fq \
-S ~/ngs_course/chipseq/results/bowtie2/H1hesc_Input_Rep1_chr12_aln_unsorted.sam

NOTE: 如果fastq文件是没有经过trim的,可以用局部比对执行soft-clipping,加上参数--local

过滤reads

首先将sam文件转为bam格式,再对bam文件进行排序,接着过滤唯一比对的reads,去除线粒体reads。 转化为bam格式 使用samtools转换格式

samtools view -h -S -b \
-o H1hesc_Input_Rep1_chr12_aln_unsorted.bam \
H1hesc_Input_Rep1_chr12_aln_unsorted.sam

对bam文件排序 对bam文件按照基因组坐标排序,可以直接使用samtools,也可以使用Sambambasambamba快速处理bam和sam文件。

sambamba sort -t 2 \
-o H1hesc_Input_Rep1_chr12_aln_sorted.bam \
H1hesc_Input_Rep1_chr12_aln_unsorted.bam 

过滤唯一比对的reads

sambamba view -h -t 2 -f bam \
-F "[XS] == null and not unmapped " \
H1hesc_Input_Rep1_chr12_aln_sorted.bam > H1hesc_Input_Rep1_chr12_aln.bam

去除PCR重复 PCR扩增和一些重复序列(如微卫星、着丝粒)会产生重复,干扰真实的富集信号,所以在call peaks前需要先去除重复,这里先用picard去除PCR重复。

java -jar picard-tools-1.119/MarkDuplicates.jar REMOVE_DUPLICATES=true I=H1hesc_Input_Rep1_chr12_aln.bam O=H1hesc_Input_Rep1_chr12_aln.dedup.bam M=H1hesc.duplicates.log

过滤线粒体reads

samtools index H1hesc_Input_Rep1_chr12_aln.dedup.bam
samtools idxstats H1hesc_Input_Rep1_chr12_aln.dedup.bam > H1hesc_Input_Rep1_chr12_aln.dedup.mitochondrial.stats
samtools view -h H1hesc_Input_Rep1_chr12_aln.dedup.bam | grep -v 'chrM' | samtools view -bS -o H1hesc.final.bam

上面给出的仅是示例代码,和参考课程不一样,实际运行需要修改相应文件。 此时就得到了唯一比对且已经去除过线粒体的比对文件,可以用于接下来的peaks calling。

参考资料:

HBC的深度NGS数据分析课程:https://github.com/hbctraining/In-depth-NGS-Data-Analysis-Course/blob/master/sessionV/lessons/02_QC_and_alignment.md

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2018-09-02

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏计算机视觉与深度学习基础

Leetcode 65 Valid Number DFA有限状态机

Validate if a given string is numeric. Some examples: "0" => true " 0.1 " =>...

27460
来自专栏大数据智能实战

tensorflow.models.rnn.rnn_cell.linear在tensorflow1.0版本之后找不到(附tensorflow1.0 API新变化)

由于版本更新关系,从原来的tensorflow低版本到升级到tensorflow1.0以上时,发现有很多API函数变化是很正常的事情,大多碰到的如: 如其中tf...

34070
来自专栏大闲人柴毛毛

张全蛋和李小花的故事——“代理模式”

举个例子: 假设张全蛋喜欢上了李小花,张全蛋准备给李小花送一些礼物,但张全蛋又非常害羞,因此买了礼物之后让王尼玛代送。王尼玛跑到李小花把礼物送给她之后就走了。 ...

404120
来自专栏点滴积累

geotrellis使用(八)矢量数据栅格化

目录 前言 栅格化处理 总结 参考链接 一、前言        首先前几天学习了一下Markdown,今天将博客园的编辑器改为Markdown,从编写博客到界面...

51970
来自专栏Jerry的SAP技术分享

Association, Composition and Aggregation in UI5, CRM, S/4HANA and C4C

比如UI5的转盘控件Carousel: 一旦转盘被析构,里面显示的page当然也没有继续存在的意义了,需要跟着被析构。

31650
来自专栏人工智能LeadAI

tensorflow读取数据-tfrecord格式

概述关于tensorflow读取数据,官网给出了三种方法: 1、供给数据:在tensorflow程序运行的每一步,让python代码来供给数据 2、从文件读取数...

1.2K60
来自专栏漏斗社区

HASH函数烧脑大作战

本期讲解一下hash函数,由于之前在比赛中做到了一题hash有关的题目,引发了此次的深(烧)度(脑)研究,本来想讲讲原理,但是太难,看得很痛苦,所以此次通过结合...

15250
来自专栏龙首琴剑庐

Java基准测试利器OpenJDK-JMH

63890
来自专栏数据结构与算法

3002 石子归并 3

 时间限制: 1 s  空间限制: 256000 KB  题目等级 : 钻石 Diamond 题解  查看运行结果 题目描述 Description 有n堆石子...

35470
来自专栏WOLFRAM

用 Wolfram 语言绘制电子轨道

16350

扫码关注云+社区

领取腾讯云代金券