MeRIP-seq 数据分析之质控、过滤、比对

原创 JunJunLab 老俊俊的生信笔记 2022-08-15

收录于合集

#100天生信/数据科学自我挑战 100 个

#MeRIP-seq 24 个

点击上方关注“公众号”

1过滤

上节讲到我们已经把数据下载下来了，总共 20 个 fastq 数据：

然后正式进入分析阶段。

创建分析环境

安装好分析需要的软件：

# 创建分析环境
$ conda create -n m6aseq python=3
# 激活环境
$ conda activate m6aseq
# 安装软件
$ conda install fastqc hisat2 samtools bedtools trim-galore cutadapt macs2 deeptools -y
# 单独安装multiqc，不然会报错后面
$ conda install -y -c conda-forge -c bioconda click multiqc

质量控制

首先用 fastqc 和 multiqc 看看数据质量和有没有 接头污染：

$ nohup mkdir -p 2.QC-data/clean-qc  && \
      fastqc -t 12 -q -o 2.QC-data/clean-qc 1.raw-data/*.gz && \
      multiqc 2.QC-data/clean-qc/* -o 2.QC-data/clean-qc/ &

所有样本的测序深度，GC 含量和重复百分比：

reads 每个碱基的质量值，看起来还是很不错的：

GC 含量的分布，理论是正态分布的，有个别样本有些异常峰：

未知碱基 N 的含量，在绿色区域内表示可以接受：

序列长度都是 150bp：

可以看到含有很多 illumina 接头的污染，我们需要去除这些接头序列：

cutadapt 过滤接头

我们使用 cutadapt 软件进行接头过滤，写个脚本批量过滤：

# 创建文件夹
$ mkdir -p 2.QC-data/trim-qc 3.tirm-data

# 写批量过滤脚本
$ vi cut.sh
#!/bin/bash
## tirm adapters
for i in SRR147656{38..47}
do
 cutadapt -a AGATCGGAAGAGC -A AGATCGGAAGAGC -m 15 -j 6 -e 0.1 -q 20 \
   -o 3.tirm-data/${i}.trimmed_1.fastq.gz \
   -p 3.tirm-data/${i}.trimmed_2.fastq.gz \
   1.raw-data/${i}_1.fastq.gz 1.raw-data/${i}_2.fastq.gz
 echo "This is sample ${i} trimming step done!"
done

# 挂后台运行
$ nohup ./cut.sh &

后台已经跑起来了：

过滤结束后，我们再质检看看质量：

$ nohup fastqc -t 12 -q -o \
        2.QC-data/trim-qc 3.tirm-data/*.gz && \
        multiqc 2.QC-data/trim-qc/* \
        -o 2.QC-data/trim-qc/ &

接头基本去除干净。

2下载索引文件

我们可以自己下载基因组文件，使用 hisat2-build 命令自己构建，也可以直接去 hisat2 官网 下载已经构建好的索引，为了方便，这里选择第二种方式，地址：http://daehwankimlab.github.io/hisat2/download/#m-musculus：

下载 UCSC 的 mm10 版本，3.54G ，我用的迅雷下的：

$ mkdir index && cd index
$ wget https://genome-idx.s3.amazonaws.com/hisat/mm10_genome.tar.gz

# 解压
$ tar -zxvf mm10_genome.tar.gz
# 查看文件
$ ls -l mm10
total 3960716
-rwxrwxrwx 1 junjun junjun  888464725 Nov 20  2015 genome.1.ht2
-rwxrwxrwx 1 junjun junjun  663195880 Nov 20  2015 genome.2.ht2
-rwxrwxrwx 1 junjun junjun       6119 Nov 20  2015 genome.3.ht2
-rwxrwxrwx 1 junjun junjun  663195875 Nov 20  2015 genome.4.ht2
-rwxrwxrwx 1 junjun junjun 1165549977 Nov 20  2015 genome.5.ht2
-rwxrwxrwx 1 junjun junjun  675339278 Nov 20  2015 genome.6.ht2
-rwxrwxrwx 1 junjun junjun          8 Nov 20  2015 genome.7.ht2
-rwxrwxrwx 1 junjun junjun          8 Nov 20  2015 genome.8.ht2
-rwxrwxrwx 1 junjun junjun       1289 Nov 20  2015 make_mm10.sh

3比对

去完接头得到我们的 clean data 后，我们可以进行比对了，使用 hisat2 比对软件，时间更快，消耗内存也较小，同时还是专门正对于转录组的比对分析软件，tophat 的升级版。

注意： 作者的测序数据是 链特异性测序，这里比对时需要加上 --rna-strandness RF 参数。

samtools 参数解释

-f 2： 提取完全匹配的序列，没有错配和缺失。
-F 256： 过滤掉多比对的序列。
-q： 提取比对质量值 30 及以上的序列。
-@： 线程数。

批量比对

比对后过滤的条件是比较严格的：

# 写个比对脚本
$ vi map.sh
#!/bin/bash
## tirm adapters
for i in SRR147656{38..47}
do
 mkdir 4.map-data
 # mapping step
 hisat2 -x index/mm10/genome -p 10 --rna-strandness RF \
 --summary-file 4.map-data/${i}.summary.txt \
 -1 3.tirm-data/${i}.trimmed_1.fastq.gz \
 -2 3.tirm-data/${i}.trimmed_2.fastq.gz \
 | samtools view -@ 10 -bS -f 2 -F 256 -q 30 \
 | samtools sort -@ 10 -o 4.map-data/${i}.high_quality_sorted.bam
done

放在后台运行：

$ nohup ./map.sh &

统计比对率

比对过程结束后,统计以下比对率,我们可以新建一个 R project 项目，后面方便管理：

创建一个 step1_summary_mapping_info.R 的脚本：

library(tidyverse)

# 设置工作路径
setwd('F:/MeRIP-seq/4.map-data')

# 文件名称
sample = list.files(pattern = '.txt')

# 提取信息
lapply(sample, function(x){
  # 读取文件
  align_info <- read.delim(x,header = F,fill = T)
  # 提取唯一比对率
  unique_map = sapply(strsplit(align_info$V1[4],split = "\\(|\\)"), '[',2)
  # 提取总比对率
  totol_map = substr(align_info$V1[15],1,6)
  # 样本名称
  name = substr(x,1,11)
  # 合并
  res = cbind(name,unique_map,totol_map)
}) %>% Reduce('rbind',.) %>%
  as.data.frame() -> mapping_info

mapping_info

          name unique_map totol_map
1  SRR14765638     76.88%    90.31%
2  SRR14765639     59.28%    69.69%
3  SRR14765640     82.85%    92.52%
4  SRR14765641     72.40%    81.10%
5  SRR14765642     82.74%    92.60%
6  SRR14765643     72.80%    81.54%
7  SRR14765644     83.06%    93.36%
8  SRR14765645     82.59%    91.09%
9  SRR14765646     79.46%    90.31%
10 SRR14765647     84.03%    91.44%

看着整体比对率还不错，SRR14765639 这个样本比对率有点低。

给 bam 文件建索引

这里顺便给 bam 文件建立索引，方便后续使用需要：

$ ls 4.map-data/*.bam | xargs -i samtools index -@ 10 {}

# 检查文件
$ ls -l 4.map-data/*.bai | cut -d ' ' -f 8,9
13:51 4.map-data/SRR14765638.high_quality_sorted.bam.bai
13:52 4.map-data/SRR14765639.high_quality_sorted.bam.bai
13:54 4.map-data/SRR14765640.high_quality_sorted.bam.bai
13:56 4.map-data/SRR14765641.high_quality_sorted.bam.bai
13:58 4.map-data/SRR14765642.high_quality_sorted.bam.bai
14:00 4.map-data/SRR14765643.high_quality_sorted.bam.bai
14:02 4.map-data/SRR14765644.high_quality_sorted.bam.bai
14:03 4.map-data/SRR14765645.high_quality_sorted.bam.bai
14:05 4.map-data/SRR14765646.high_quality_sorted.bam.bai
14:07 4.map-data/SRR14765647.high_quality_sorted.bam.bai