生信 | scATAC-Seq数据获取（超快sra转fq）

2021-07-12 本文已影响0人生信卷王

部分代码参考自10x genomics官网，但绝大多数代码融入了我自己的思考

1、上游分析（软件安装 + 数据获取 + call peak）

由于10X封装好了流程，所以软件安装和call peak就很简单了，因此重要的就是如何获取数据？、如何超快速从sra跑出来R1234？，不过还是解释一下整体流程吧，不想看的可以直接往下翻！

1.1 软件安装

1）服务器硬件与软件要求

硬件要求：>8核CPU，>64G运行内存，>1TB存储空间
软件要求：bcl2fastq版本 ≥ 2.20【官网链接(要注册)】【从conda安装(推荐)】

# 先激活任意一个 conda 环境，然后再执行下面的代码
conda install -c freenome bcl2fastq -y
# 检查是否安装成功
bcl2fastq -h

2） Cell Ranger ATAC - 2.0.0软件安装（1.7GB）

因为经常更新，所以去【官网】看一看是不是2.0.0版本了，至少在我更新的时候还是的

wget -O cellranger-atac-2.0.0.tar.gz "https://cf.10xgenomics.com/releases/cell-atac/cellranger-atac-2.0.0.tar.gz?Expires=1626103123&Policy=eyJTdGF0ZW1lbnQiOlt7IlJlc291cmNlIjoiaHR0cHM6Ly9jZi4xMHhnZW5vbWljcy5jb20vcmVsZWFzZXMvY2VsbC1hdGFjL2NlbGxyYW5nZXItYXRhYy0yLjAuMC50YXIuZ3oiLCJDb25kaXRpb24iOnsiRGF0ZUxlc3NUaGFuIjp7IkFXUzpFcG9jaFRpbWUiOjE2MjYxMDMxMjN9fX1dfQ__&Signature=VODjodUMuNGP51unCTQHn6ONd8dI~tepKSW0qDcxv2n2LKQoCMxKoTD6Vdkidwl7F~a0XAnc1~ji5V4wgkQgeY5~gbqf~ZMNpAxf3AXU9XreMwsP7izc-EDOelCbcUkK8n3ynfzXSKvq9qR008LjPwIGUAN0prS9xJNFkPD0Dq7GiK~tKZ5-g~20YEwpx--OMeIpzEvVlq-xpKNVy8qjulshKgslmcOPs7lHxv9bI04Xd8XFS5CE9DoKKD~2bCyXTXWYSRAhWLsXDdozGzVaSVJwh8sbarQzRaot79GCS~xQPkmZJmnXuKgHFiwUDmFenOtpkYh2gqpo1nEfaFBlTQ__&Key-Pair-Id=APKAI7S6A5RYOXBWRPDA"
tar -xzvf cellranger-atac-2.0.0.tar.gz
# 永久添加到环境变量
echo 'export PATH=/你的安装路径/cellranger-atac-2.0.0:$PATH' >> ~/.bashrc
source ~/.bashrc
# 检查是否安装成功，需要 3分钟，最后提示Pipestance completed successfully!
cellranger-atac testrun --id=tiny
rm cellranger-atac-2.0.0.tar.gz

1.2 数据获取

如果你已经有了I1,R1,R2,R3以下内容可以跳过，主要阐述分析的I1,R1,R2,R3数据是怎么得到的。来源有两个： $\color{red}{①\ \ 下机的BCL转fastq}$ $\color{red}{②\ \ 从SRA/ENA上下载fastq}$

1）数据来源一：对下机数据base calling files（BCL）转为fastq文件（做生信的跳过）
软件：cellranger mkfastq ：它借鉴了Illumina的bcl2fastq ，可以将一个或多个lane中的混样测序样本按照index标签生成样本对应的fastq文件。即对下机数据base calling files转为fastq文件。

下载示例数据（233MB）

# 下载Illumina® BCL的测序文件，大小为233MB
wget https://cf.10xgenomics.com/supp/cell-atac/cellranger-atac-tiny-bcl-1.0.0.tar.gz
# 下载样本注释文件，也可以自己手动创造，见下文格式
wget https://cf.10xgenomics.com/supp/cell-atac/cellranger-atac-tiny-bcl-simple-1.0.0.csv
tar -zvxf cellranger-atac-tiny-bcl-1.0.0.tar.gz
rm cellranger-atac-tiny-bcl-1.0.0.tar.gz

可以看一下样本注释文件的格式

head cellranger-atac-tiny-bcl-simple-1.0.0.csv
# 三列，逗号分割
Lane,Sample,Index
1,test_sample,SI-NA-C1

运行程序——mkfastq

cellranger-atac mkfastq --id=tiny-bcl \
                     --run=cellranger-atac-tiny-bcl-1.0.0 \
                     --csv=cellranger-atac-tiny-bcl-simple-1.0.0.csv \
                     --delete-undetermined

--id：输出文件夹的名字
--run：测序文件夹的名字
--csv：样本注释文件
--delete-undetermined：删除输出文件中的undetermined文件

查看运行结果

ls -lh tiny-bcl/outs/fastq_path/HJN3KBCX2/test_sample

一个样本生成四个文件，分别为I1, R1,R2,R3，依次代表着sample index,read1,barcode,read2
大多数情况下，我们直接会拿到一个样本的四个文件I1, R1,R2,R3或者其中的R1,R3，所以上面使用mkfastq对Illumina®BCL的测序文件进行的处理可以跳过，直接进行下面的分析即可

$\color{red}{2）数据来源二：从SRA/ENA上下载（生信砖家请重点关注）}$

这里面的门道感觉还比较多，经过我3天的探索，基本上摸清楚了如何快速得到scATAC-Seq产生的I1, R1,R2,R3的四个文件，注意是快速！请听我慢慢道来（又在废话）
软件：【Aspera Connect】 + 【SRA Toolkit】【帮助文档】
以前我是从来不使用SRA Toolkit，但是现在想获取scATAC产生的数据不得不用到，没办法，人在屋檐下哪能不低头？另外我还想说一点的是SRA Toolkit目前依旧支持 ascp 下载，不知道以前听谁说不可以的
软件安装（均不推荐使用conda安装）

wget https://d3gcli72yxqn2z.cloudfront.net/connect_latest/v4/bin/ibm-aspera-connect-3.11.2.63-linux-g2.12-64.tar.gz
tar -zvxf ibm-aspera-connect-3.11.2.63-linux-g2.12-64.tar.gz
sh ibm-aspera-connect-3.11.2.63-linux-g2.12-64.sh
# 永久添加到环境变量
echo 'export PATH=~/.aspera/connect/bin:$PATH' >> ~/.bashrc

wget https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.9.6/sratoolkit.2.9.6-ubuntu64.tar.gz
tar -zvxf sratoolkit.2.9.6-ubuntu64.tar.gz
echo 'export PATH=/你的安装路径/sratoolkit.2.9.6-ubuntu64/bin:$PATH' >> ~/.bashrc

source ~/.bashrc

下载数据
首先不同于传统的NGS测序数据，如果你直接下载_1.fastq.gz和_2.fastq.gz由于缺少barcode文件，cellranger ATAC是无法运行的，所以只能使用SRA Toolkit将原始sra文件下载下来，然后再分解为单独的文件。
好！重点来了！既然要分解那就要选择块一点的方式，传统的fastq-dump我是不敢再用了，巨慢。因此有没有替代品呢？那是当然！同样来自SRA Toolkit，叫做fasterq-dump。顾名思义，比fastq-dumpfaster的分解软件，下面正式开始
以下载SRR13450125.sra为例

prefetch SRR13450125 -O ./

-O：大写的O，表示存放在指定的文件夹下
使用prefetch是极度看脸的，因为你设置的所有想让他使用ascp下载的参数后，他依旧可能使用https下载，因此佛系一点吧，https也不算太慢。

使用fasterq-dump分解

fasterq-dump -p --include-technical -S -e 10 -O ./ SRR13450125.sra

-p：显示进程，非常推荐
--include-technical：这是关键点，如果不加这个参数是不会分解出来index和barcode文件的，这就不同于fastq-dump的--split-files
-S：将分解出来的序列存放到不同的文件中，如I1, R1,R2,R3四个文件
-e：设置线程数
-O：大写的O，表示存放在指定的文件夹下

fasterq-dump是真的快啊，不到2分钟，就分解完了，不过不足之处就是输出不是压缩文件，毕竟压缩也是需要时间的，因此用内存换速度，我觉的挺值的

修改fastq的名字
因为CellRanger ATAC只能识别他pipeline规定的名字格式，因此需要改名，格式如下

[Sample Name]_S1_L00[Lane Number]_[Read Type]_001.fastq
--------------------------------------------------------------
[Sample Name]：自己对样本的命名
[Lane Number]：同一个样本不同测序泳道的序号
[Read Type]：四种类型，如下
I1: Sample index read (这个文件可要可不要)
R1: Read 1
R2: barcode（这个文件必须要有）
R3: Read 2
--------------------------------------------------------------
举例如下：
SRR13450125_S1_L001_I1_001.fastq
SRR13450125_S1_L001_R1_001.fastq
SRR13450125_S1_L001_R2_001.fastq
SRR13450125_S1_L001_R3_001.fastq

可是，我们怎么知道上面的1,2,3,4怎么对应I1, R1,R2,R3啊？方法很简单，在浏览器中输入下面的网址查询一下文件的基本信息，需要查询其他的修改run=SRR13450125即可

https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR13450125

熟悉barcode的应该知道它为16bp，因此其中的L=16就代表着barcode文件，也就是R2，所以SRR13450125.sra_3.fastq就可以改名为SRR13450125_S1_L001_R2_001.fastq。而R1和R2命名给SRR13450125.sra_2.fastq或者SRR13450125.sra_4.fastq都可以，不放心的可以head一些文件看一看，这里我跑了一下fastqc进一步验证了我的观点

所以根据以上规则改名：

mv SRR13450125.sra_1.fastq SRR13450125_S1_L001_I1_001.fastq
mv SRR13450125.sra_2.fastq SRR13450125_S1_L001_R1_001.fastq
mv SRR13450125.sra_3.fastq SRR13450125_S1_L001_R2_001.fastq
mv SRR13450125.sra_4.fastq SRR13450125_S1_L001_R3_001.fastq

作者不想让你分析系列
如果你探索的足够深入的话，你还会发现一些奇怪的东西，比如SRR14539571

这种就直接放弃吧，我看了一下他是把R1(L=51)R2(L=16)R3(L=51)序列合并了，很难搞，感兴趣的同学可以自己把barcode提取出来继续分析，这里不再展开

$ head SRR14539571.sra_2.fastq
@SRR14539571.sra.1 NB501658:300:HVW2LBGXF:1:11101:15875:1034 length=118
GTCCANCCATTTGTCTAAGAATTTCATGAATTCGTTGTTTCTAATAGCTAATGGCCGATGCGATGGACTCTTCCCATTGATGCCCGCNTAGGTNATNCTCTGCTACATATGCATCTAN
+SRR14539571.sra.1 NB501658:300:HVW2LBGXF:1:11101:15875:1034 length=118
AAAA/#AEE6/E6EE/EAEAEEAEEEEEEEEE6/EEEE/EEEEE/EAEAA/AAAAAEEAAEEEAEA/6AAAA/A/A6EEEE///E/E#EEEEE#AE#/E/A///EAEEEE//EA<AE#
@SRR14539571.sra.2 NB501658:300:HVW2LBGXF:1:11101:10442:1035 length=118
GTTTANTGGCTAGCTCATCCAGACCCCTCAACTAAACTAGTTTTCTATTGGCCCTTGCCTACAGGCATCGCTATTCTCGGCCCTTCCNGTTCAACANCTAGCCCTGGGTTCCCGAAGG
+SRR14539571.sra.2 NB501658:300:HVW2LBGXF:1:11101:10442:1035 length=118
AA6AA#EEE/E//E/EAEEE/E/AEEEEEE6EEEE6EE6EEEEEE/EEEEEA/A/AEAE/AEAEEAAAAAAAEAEE///EEEEE//A#A/E/EEE/#AEE/E6/6E//<A<EAAE/AE

一个样本多个Lane系列
如果还要继续探索的，还存在另外一种情况，就是一个样本有多个Lane，因此分析的时候要把原始文件下载全了再分析，就比如说上面下载的SRR13450125，他其实还有个孪生兄弟SRR13450126，只不过测序的时候泳道不同，所以下载的时候一定要下载全，虽然不全也能继续分析，只不过就不能复原文献内容了。这种情况在传统ATAC-Seq中也存在，这里不再展开叙述了

终于把下载数据折腾完了

1.3 下载【比对索引】或者【自己构建索引】

这里就不得不吐槽10X了，你整那么大的索引干什么，下载起来贼费劲。

下载适配于Cell Ranger ATAC小鼠的参考基因组（13GB，解压后19GB）。没办法，就是这么大，但下载总比自己构建的快且用的安心吧

wget https://cf.10xgenomics.com/supp/cell-atac/refdata-cellranger-arc-mm10-2020-A-2.0.0.tar.gz
tar -zvxf refdata-cellranger-arc-mm10-2020-A-2.0.0.tar.gz
rm refdata-cellranger-arc-mm10-2020-A-2.0.0.tar.gz #快点删了

番外篇：自己构建适配于Cell Ranger ATAC的参考基因组——mkref
目前官网上提供的参考基因组只有GRCh38与mm10，因此做其他动物或植物的需要自己构建索引。但必须同时具有gtf和genome文件，所以不想下载或官网上没研究的物种的，也可以自己构建

# 首先准备一个配置文件
vi config
# 在配置文件中输入以下内容后保存（字典格式）
{
    organism: "mouse"
    genome: ["GRCh38"]
    input_fasta: ["/public/bioinfoYu/genome/GRCm38.assembly.genome.fa"]
    input_gtf: ["/public/bioinfoYu/genome/GRCm38.annotation.gtf"]
    non_nuclear_contigs: ["chrM"]
    input_motifs: "/public/bioinfoYu/genome/motifs.pfm"
}
# 然后执行程序
cellranger-atac mkref --config=config

organism：可选参数，指定你想定义的物种名
genome：必选参数，输出文件所存放的文件夹名称
input_fasta：必选参数，物种的基因组文件
input_gtf：必选参数，物种的gtf注释文件
non_nuclear_contigs：构建索引的时候忽略的染色体，chrM指线粒体上的染色体，指定多个用逗号分割，如：["chrM","chrC"]
input_motifs:：可选参数，指定jaspar格式的motif文件，推荐加上啊，可以去扣扣群559758504下载hg38或者mm10的motifs文件

1.4 call peak

这里就以我们下载的SRR13450125和SRR13450126为例吧，因为二者是孪生兄弟，因此分析的时候一起分析，所以改名的时候不能说改为SRR13450125_S1...和SRR13450126_S1...，要命名成同样的前缀，这里我统一命名为armstrong_IGO_09847_1，和作者保持一致，如下：

开始分析，cellranger-atac count的原理就是单纯的将两个lane文件合并后分析，后面我们看一下差别

cellranger-atac count --id=armstrong_result \
                        --reference=refdata-cellranger-arc-mm10-2020-A-2.0.0 \
                        --fastqs=mm10/ \
                        --sample=armstrong_IGO_09847_1 \
                        --localcores=8 \
                        --localmem=64

--id：输出的文件夹名称
--reference：下载或自己构建的参考基因组索引文件夹
--fastqs：存放样本的文件夹。路径前面一定不要有【/】，否则报错
--sample：样本名，也就是_S1前面的内容，我们前面改名为了armstrong_IGO_09847_1_S1...因此这个地方就是armstrong_IGO_09847_1
--lanes：指定分析一个样本下的单个lane，如--lanes=1则只分析L001的内容。如果不使用，则默认分析一个样本下的所有lanes，这里我们不使用，则同时分析L001和L002
--localcores：默认占满服务器的所有核CPU，因此推荐设置的小一点，如8核
--localmem：请调整参数限制占用的运行内存

注意：想了解更多--fastqs，--sample，--lanes如何配合使用，请参考【官网】或【中文解释】

为了比较差异，我单独跑了一下L001和L002的call peak结果，差异比较放到结果解读后面吧

for i in 1 2
do
cellranger-atac count --id=armstrong_result_${i} \
                        --reference=refdata-cellranger-arc-mm10-2020-A-2.0.0 \
                        --fastqs=mm10/test \
                        --sample=armstrong_IGO_09847_1 \
                        --lanes=${i}
                        --localcores=8 \
                        --localmem=64
done

1.5 结果解读

Outputs:
- Per-barcode fragment counts & metrics:        /home/bioinfoYu/armstrong_result/outs/singlecell.csv
- Position sorted BAM file:                     /home/bioinfoYu/armstrong_result/outs/possorted_bam.bam
- Position sorted BAM index:                    /home/bioinfoYu/armstrong_result/outs/possorted_bam.bam.bai
- Summary of all data metrics:                  /home/bioinfoYu/armstrong_result/outs/summary.json
- HTML file summarizing data & analysis:        /home/bioinfoYu/armstrong_result/outs/web_summary.html
- Bed file of all called peak locations:        /home/bioinfoYu/armstrong_result/outs/peaks.bed
- Raw peak barcode matrix in hdf5 format:       /home/bioinfoYu/armstrong_result/outs/raw_peak_bc_matrix.h5
- Raw peak barcode matrix in mex format:        /home/bioinfoYu/armstrong_result/outs/raw_peak_bc_matrix
- Directory of analysis files:                  /home/bioinfoYu/armstrong_result/outs/analysis
- Filtered peak barcode matrix in hdf5 format:  /home/bioinfoYu/armstrong_result/outs/filtered_peak_bc_matrix.h5
- Filtered peak barcode matrix in mex format:   /home/bioinfoYu/armstrong_result/outs/filtered_peak_bc_matrix
- Barcoded and aligned fragment file:           /home/bioinfoYu/armstrong_result/outs/fragments.tsv.gz
- Fragment file index:                          /home/bioinfoYu/armstrong_result/outs/fragments.tsv.gz.tbi
- Filtered tf barcode matrix in hdf5 format:    /home/bioinfoYu/armstrong_result/outs/filtered_tf_bc_matrix.h5
- Filtered tf barcode matrix in mex format:     /home/bioinfoYu/armstrong_result/outs/filtered_tf_bc_matrix
- Loupe Browser input file:                     /home/bioinfoYu/armstrong_result/outs/cloupe.cloupe
- csv summarizing important metrics and values: /home/bioinfoYu/armstrong_result/outs/summary.csv
- Annotation of peaks with genes:               /home/bioinfoYu/armstrong_result/outs/peak_annotation.tsv
- Peak-motif associations:                      /home/bioinfoYu/armstrong_result/outs/peak_motif_mapping.bed
 
Pipestance completed successfully!