miRNA生物信息数据分析流程初探(二)

最后发布时间:2022-11-11 16:41:02 浏览量:

上篇文章中我们描述了miRNA的基本分析流程,本文中我们将对miRNA分析流程中每一个步骤的输出文件进行解释。

软件安装

在snakemake流程中可以使用conda去管理需要的软件,详细的命令使用见: Running Snakemake in conda env, 下面介绍其基本用法:

查看snakemakeconda的环境

snakemake  --directory .test   --list-conda-envs

结果显示如下

environment     container       location
../workflow/envs/env.yaml               .snakemake/conda/369d9fb0f93ac536da2fc3c2969c779d

使用conda的环境

snakemake -p --directory .test resources/mirna.fasta   --use-conda --conda-frontend mamba -c5 

结果显示如下

Building DAG of jobs...
Creating conda environment ../workflow/envs/env.yaml...
Downloading and installing remote packages.

在第一次使用conda环境的时候,会首先安装,该过程由于国内网速的原因,比较耗时间。这里我们可以手动创建conda的环境,具体的命令如下:

mamba env create  --file  workflow/envs/env.yaml --prefix .snakemake/conda/369d9fb0f93ac536da2fc3c2969c779d

注意--prefix后面的内容是--list-conda-envs中查找到的,并且workflow/envs/env.yaml文件名称及内容一旦修改,369d9fb0f93ac536da2fc3c2969c779d这个hash值就会更改

如果你此时运行snakemake命令,仍旧会创建环境,此时需要运行下面命令创建一个空的文件夹

touch .test/.snakemake/conda/369d9fb0f93ac536da2fc3c2969c779d/env_setup_end

数据库下载

miRNA数据下载

miRNA的数据库地址: https://www.mirbase.org/ftp.shtml
可以使用下面命令下载

 snakemake  -p --directory .test resources/mirna.fasta   --use-conda --conda-frontend mamba -c5

此时会得到三个文件

.
├── mirna.fasta
├── mirna_hairpin.fasta
└── mirna_mature.fasta

#  mirna_hairpin.fasta
>hsa-let-7a-1 MI0000060 Homo sapiens let-7a-1 stem-loop
TGGGATGAGGTAGTAGGTTGTATAGTTTTAGGGTCACACCCACCACTGGGAGATAACTAT
ACAATCTACTGTCTTTCCTA
>hsa-let-7a-2 MI0000061 Homo sapiens let-7a-2 stem-loop
AGGTTGAGGTAGTAGGTTGTATAGTTTAGAATTACATCAAGGGAGATAACTGTACAGCCT
CCTAGCTTTCCT
>hsa-let-7a-3 MI0000062 Homo sapiens let-7a-3 stem-loop
GGGTGAGGTAGTAGGTTGTATAGTTTGGGGCTCTGCCCTGCTATGGGATAACTATACAAT
CTACTGTCTTTCCT
>hsa-let-7b MI0000063 Homo sapiens let-7b stem-loop

# head mirna.fasta
>hsa-let-7a-5p MIMAT0000062 Homo sapiens let-7a-5p
TGAGGTAGTAGGTTGTATAGTT
>hsa-let-7a-3p MIMAT0004481 Homo sapiens let-7a-3p
CTATACAATCTACTGTCTTTC
>hsa-let-7a-2-3p MIMAT0010195 Homo sapiens let-7a-2-3p
CTGTACAGCCTCCTAGCTTTCC
>hsa-let-7b-5p MIMAT0000063 Homo sapiens let-7b-5p
TGAGGTAGTAGGTTGTGTGGTT
>hsa-let-7b-3p MIMAT0004482 Homo sapiens let-7b-3p
CTATACAACCTACTGCCTTCCC

# mirna_mature.fasta 是上面两个文件合并后的文件
miRBase Sequence

hairpin.fa: Fasta format sequences of all miRNA hairpins(发夹结构)
mature.fa: Fasta format sequences of all mature(成熟) miRNA sequences

生信小木屋

图a: 前体(precursor)miRNA序列发夹结构(hairpin conformation)的示意图,成熟的miRNA位于茎(蓝色突出显示)及其相应的种子区域(红色)SL:mature_microrna
图b: 典型的miRNA–mRNA相互作用,其中每个结合位点通过6至7个Watson–Crick配对(垂直线)与成熟miRNA的种子区(第2至第8个5′核苷酸)配对

piRNA数据下载

piRNA数据下载地址SL:pirna_base: http://bigdata.ibp.ac.cn/piRBase/download/v3.0/fasta/hsa.v3.0.fa.gz
可以使用下面命令下载

snakemake  -p --directory .test resources/pirna.fasta   --use-conda --conda-frontend mamba -c5

此时会得到一个fasta文件

.
└── pirna.fasta

# head pirna.fasta
>piR-hsa-1
AGAAGACATTCGTGGAGGCGTC
>piR-hsa-2
ACGCCTCCACGTAGTGTCTT
>piR-hsa-3
ACGCCTCCACGAGTGTCTT
>piR-hsa-4
ACCAAGGACTTCATCTTCTCAGAGCTGCTGT
>piR-hsa-5
TCATCTTCTCAGAGCTGCTGTCCAACCTGC

建立index

miRNA建立index

snakemake  --directory .test  resources/mirna.1.bt2  --use-conda --conda-frontend mamba -c5   -p

得到以下文件

.
├── logs
├── mirna.1.bt2
├── mirna.2.bt2
├── mirna.3.bt2
├── mirna.4.bt2
├── mirna.rev.1.bt2
└── mirna.rev.2.bt2

piRNA建立index

使用下面命令建立piRNA的index

snakemake  --directory .test  resources/pirna.1.bt2  --use-conda --conda-frontend mamba -c5   -p

得到以下文件

.
├── pirna.1.bt2
├── pirna.2.bt2
├── pirna.3.bt2
├── pirna.4.bt2
├── pirna.rev.1.bt2
└── pirna.rev.2.bt2

下机fastq数据质控

merge_fastq_files

合并单个样本所有的fastq文件到一个文件中

snakemake  --directory .test  results/sampleA/fastq/sampleA_untrimmed.fastq.gz  --use-conda --conda-frontend mamba -c5   -p

合并所有sampleAfast1文件得到文件sampleA_untrimmed.fastq.gz

# zcat sampleA_untrimmed.fastq.gz | head

@Simulated1
CGAGGAACTGTAGGCACCATCAATAAAACCCTAGAC [AGATCGGAAGAGCACACGTCTGAAC] TCCAGTCACTTAGG
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEAEEEEEEEAEEEEAEAEEAEEEEEEEEEEEE
@Simulated2
CAGTTTTGTCAGGTAAACTGTAGGCACCATCAATCCCAGTCGCTTC [AGATCGGAAGAGCACACGTCTGAAC] TCCA
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEAEE<EAAAEEEEEEEA/EEEEEEAEEEEEAE
@Simulated3
TGCTGCTTTATTGTTATTGCTGCTTGATTGTTAACTGTAGGCACCATCAATATGAAGACGAATAGATCGGAAGAG

上面括号中的序列是Illumina的通用接头序列

trim_illumina_adapters

从reads的3'端去除Illumina的通用接头序列(Adapter sequences)
通用接头序列是: AGATCGGAAGAGCACACGTCTGAAC

snakemake  --directory .test  results/sampleA/fastq/sampleA_illumina_trimmed.fastq.gz  --use-conda --conda-frontend mamba -c5   -p

结果得到文件sampleA_illumina_trimmed.fastq.gz

# zcat sampleA_illumina_trimmed.fastq.gz | head

@Simulated1
CGAGGAACTGTAGGCACCATCAAT [AAAACCCTAGAC]
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEE
@Simulated2
CAGTTTTGTCAGGTAAACTGTAGGCACCATCAAT [CCCAGTCGCTTC]
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEAE
@Simulated3
TGCTGCTTTATTGTTATTGCTGCTTGATTGTTAACTGTAGGCACCATCAAT [ATGAAGACGAAT]

上面括号中的序列为UMI

trim_umi_sequences

从reads的3'端去除UMI(Unique Molecular Identifier)序列,并将其添加到reads的头部
其中UMI的pattern为: NNNNNNNNNNNN

snakemake  --directory .test  results/sampleA/fastq/sampleA_umi_trimmed.fastq.gz  --use-conda --conda-frontend mamba -c5   -p 

结果得到文件sampleA_umi_trimmed.fastq.gz

# zcat sampleA_umi_trimmed.fastq.gz | head

@Simulated1_AAAACCCTAGAC
CGAGGAACTGTAGGCACCATCAAT
+
AAAAAEEEEEEEEEEEEEEEEEEE
@Simulated2_CCCAGTCGCTTC
CAGTTTTGTCAGGTAAACTGTAGGCACCATCAAT
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
@Simulated3_ATGAAGACGAAT
TGCTGCTTTATTGTTATTGCTGCTTGATTGTTAACTGTAGGCACCATCAAT

关于UMI在3'端的的原因,可以看下图,在illumina测序中reads是从5'向3'端测的。UMI是在adapter之后的,因此去除adapter后紧跟着的12bp为UMI的序列

生信小木屋

可以将UMI集成到Y-adapter连接中,而无需中断标准工作流程。在文库准备期间,原始文库准备套件的adapter简单地由UMI adadpter代替SL:umi

trim_5_adapter

从reads 5'端去除接头

snakemake  --directory .test  results/sampleA/fastq/sampleA_5adapter_trimmed.fastq.gz  --use-conda --conda-frontend mamba -c5   -p

结果得到文件sampleA_5adapter_trimmed.fastq.gz

# zcat sampleA_5adapter_trimmed.fastq.gz | head

@Simulated1_AAAACCCTAGAC
CGAGGAACTGTAGGCACCATCAAT
+
AAAAAEEEEEEEEEEEEEEEEEEE
@Simulated2_CCCAGTCGCTTC
CAGTTTTGTCAGGTAAACTGTAGGCACCATCAAT
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
@Simulated3_ATGAAGACGAAT
TGCTGCTTTATTGTTATTGCTGCTTGATTGTTAACTGTAGGCACCATCAAT

trim_3_adapter

从reads 3'端去除接头,丢弃未修剪的reads并执行一些质量调整

snakemake  --directory .test  results/sampleA/fastq/sampleA_trimmed.fastq.gz  --use-conda --conda-frontend mamba -c5   -p

结果得到文件sampleA_trimmed.fastq.gz

# zcat sampleA_trimmed.fastq.gz | head

@Simulated3_ATGAAGACGAAT
TGCTGCTTTATTGTTATTGCTGCTTGATTGTT
+
AAAAAEEEEEEEEEEEEAEEEEEEEEAEEAEE
@Simulated19_ATTCCCGCCGCG
GTTGATTCTGTGGGTGGTGGTGCAT
+
/AAAAAEEEEE/EEEEEAEEEEEEE
@Simulated31_TGACCCATCCCC
GATTTATGAAAGACGAACTACTGCCAAAG

比对reads到数据库

map_mirna_sequences

将reads比对到mirBase数据库序列,直接转换为排序的.bam文件,并生成.bai的索引文件。没有比对上的reads将保存为新的fastq 文件,随后与piRNA的数据库比对。

snakemake  --directory .test  results/sampleA/bam/sampleA_mirna.bam  --use-conda --conda-frontend mamba -c5   -p

得到3个结果文件

sampleA_mirna.bai
sampleA_mirna.bam
sampleA_mirna_unmapped.fastq.gz

map_pirna_sequences

将reads比对到piRNAdb数据库序列,直接转换为排序的.bam文件,并生成.bai的索引文件。有比对上的reads将保存为新的fastq 文件,随后与参考基因组比对(本文没有执行该步骤)。
这里使用的序列是miRNA没有mapping上的文件:sampleA_mirna_unmapped.fastq.gz

snakemake  --directory .test  results/sampleA/bam/sampleA_pirna.bam  --use-conda --conda-frontend mamba -c5   -p

得到3个结果文件

sampleA_pirna.bai
sampleA_pirna.bam
sampleA_pirna_unmapped.fastq.gz

合并bam文件

snakemake  --directory .test  results/sampleA/bam/sampleA_merged.bam  --use-conda --conda-frontend mamba -c5   -p

得到两个文件

sampleA_merged.bai
sampleA_merged.bam

deduplicate_reads

使用UMItools dedup命令根据先前识别的UMI序列,去除重复的reads

snakemake  --directory .test  results/sampleA/bam/sampleA_merged_dedup.bam  --use-conda --conda-frontend mamba -c5   -p

得到两个文件

sampleA_merged_dedup.bai
sampleA_merged_dedup.bam

获得count表

snakemake  --directory .test  results/sampleA/sampleA_counts_dedup.tsv  --use-conda --conda-frontend mamba -c5   -p

得到文件sampleA_counts_dedup.tsv

# head sampleA_counts_dedup.tsv

Molecule        Count
hsa-let-7f-5p   1
hsa-miR-16-5p   5
hsa-miR-21-5p   1
hsa-miR-24-3p   1
hsa-miR-28-3p   1
hsa-miR-92a-3p  1
hsa-miR-93-5p   3
hsa-miR-221-3p  1
hsa-miR-223-3p  1

运行所有样本

snakemake  --directory .test    --use-conda --conda-frontend mamba -c5   -p

最终得到两张表,其中counts.tsv没有运行deduplicate_readssampleA_counts_dedup两条命令。

# head counts.tsv
Molecule        sampleA sampleB sampleC
hsa-let-7a-1    0.0     0.0     4.0
hsa-let-7a-2    4.0     4.0     0.0
hsa-let-7a-3    4.0     0.0     0.0
hsa-let-7b      12.0    4.0     4.0
hsa-let-7b-5p   0.0     4.0     0.0
hsa-let-7f-1    0.0     0.0     4.0
hsa-let-7f-5p   4.0     0.0     0.0
hsa-miR-103a-3p 0.0     0.0     4.0
hsa-miR-125a-5p 4.0     0.0     0.0

# counts_deduplicated.tsv
Molecule        sampleA sampleB sampleC
hsa-let-7a-1    0.0     0.0     1.0
hsa-let-7a-2    1.0     1.0     0.0
hsa-let-7a-3    1.0     0.0     0.0
hsa-let-7b      3.0     1.0     1.0
hsa-let-7b-5p   0.0     1.0     0.0
hsa-let-7f-1    0.0     0.0     1.0
hsa-let-7f-5p   1.0     0.0     0.0
hsa-miR-103a-3p 0.0     0.0     1.0
hsa-miR-125a-5p 1.0     0.0     0.0