部分代碼參考自10x genomics官網(wǎng)尽纽，但絕大多數(shù)代碼融入了我自己的思考

1夏志、上游分析（軟件安裝 + 數(shù)據(jù)獲取 + call peak）

• 由于10X封裝好了流程乃坤，所以軟件安裝和call peak就很簡單了，因此重要的就是如何獲取數(shù)據(jù)沟蔑？湿诊、如何超快速從sra跑出來R1234？瘦材，不過還是解釋一下整體流程吧厅须，不想看的可以直接往下翻！

1.1 軟件安裝

1） 服務(wù)器硬件與軟件要求

• 硬件要求：>8核CPU食棕，>64G運行內(nèi)存朗和，>1TB存儲空間
• 軟件要求：bcl2fastq版本 ≥ 2.20【官網(wǎng)鏈接(要注冊)】【從conda安裝(推薦)】

# 先激活任意一個 conda 環(huán)境，然后再執(zhí)行下面的代碼
conda install -c freenome bcl2fastq -y
# 檢查是否安裝成功
bcl2fastq -h

2） Cell Ranger ATAC - 2.0.0軟件安裝（1.7GB）

• 因為經(jīng)常更新簿晓，所以去【官網(wǎng)】看一看是不是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
# 永久添加到環(huán)境變量
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 數(shù)據(jù)獲取

• 如果你已經(jīng)有了I1,R1,R2,R3以下內(nèi)容可以跳過镀层，主要闡述分析的I1,R1,R2,R3數(shù)據(jù)是怎么得到的。來源有兩個：

1）數(shù)據(jù)來源一：對下機數(shù)據(jù)base calling files（BCL）轉(zhuǎn)為fastq文件（做生信的跳過）
軟件：cellranger mkfastq ： 它借鑒了Illumina的bcl2fastq 皿曲，可以將一個或多個lane中的混樣測序樣本按照index標簽生成樣本對應(yīng)的fastq文件唱逢。即對下機數(shù)據(jù)base calling files轉(zhuǎn)為fastq文件。

• 下載示例數(shù)據(jù)（233MB）

# 下載Illumina? BCL的測序文件屋休，大小為233MB
wget https://cf.10xgenomics.com/supp/cell-atac/cellranger-atac-tiny-bcl-1.0.0.tar.gz
# 下載樣本注釋文件坞古，也可以自己手動創(chuàng)造，見下文格式
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文件

• 查看運行結(jié)果

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

• 一個樣本生成四個文件痪枫，分別為I1, R1,R2,R3织堂，依次代表著sample index,read1,barcode,read2
  大多數(shù)情況下，我們直接會拿到一個樣本的四個文件I1, R1,R2,R3或者其中的R1,R3奶陈，所以上面使用mkfastq對Illumina?BCL的測序文件進行的處理可以跳過易阳，直接進行下面的分析即可

• 這里面的門道感覺還比較多，經(jīng)過我3天的探索吃粒，基本上摸清楚了如何快速得到scATAC-Seq產(chǎn)生的I1, R1,R2,R3的四個文件潦俺，注意是快速！請聽我慢慢道來（又在廢話）

• 軟件： 【Aspera Connect】 + 【SRA Toolkit】 【幫助文檔】
  以前我是從來不使用SRA Toolkit徐勃，但是現(xiàn)在想獲取scATAC產(chǎn)生的數(shù)據(jù)不得不用到事示，沒辦法，人在屋檐下哪能不低頭僻肖？另外我還想說一點的是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
# 永久添加到環(huán)境變量
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

• 下載數(shù)據(jù)
  首先不同于傳統(tǒng)的NGS測序數(shù)據(jù)，如果你直接下載_1.fastq.gz和_2.fastq.gz由于缺少barcode文件臀脏，cellranger ATAC是無法運行的劝堪，所以只能使用SRA Toolkit將原始sra文件下載下來，然后再分解為單獨的文件谁榜。
  好幅聘！重點來了！既然要分解那就要選擇塊一點的方式窃植，傳統(tǒng)的fastq-dump我是不敢再用了帝蒿，巨慢。因此有沒有替代品呢巷怜？那是當然葛超！同樣來自SRA Toolkit，叫做fasterq-dump延塑。顧名思義绣张，比fastq-dumpfaster的分解軟件，下面正式開始

• 以下載SRR13450125.sra為例

prefetch SRR13450125 -O ./

-O：大寫的O关带，表示存放在指定的文件夾下
使用prefetch是極度看臉的侥涵，因為你設(shè)置的所有想讓他使用ascp下載的參數(shù)后，他依舊可能使用https下載宋雏，因此佛系一點吧芜飘，https也不算太慢。

• 使用fasterq-dump分解

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

-p：顯示進程磨总，非常推薦
--include-technical：這是關(guān)鍵點嗦明，如果不加這個參數(shù)是不會分解出來index和barcode文件的，這就不同于fastq-dump的--split-files
-S：將分解出來的序列存放到不同的文件中蚪燕，如I1, R1,R2,R3四個文件
-e：設(shè)置線程數(shù)
-O：大寫的O娶牌，表示存放在指定的文件夾下

• fasterq-dump是真的快啊奔浅，不到2分鐘，就分解完了诗良，不過不足之處就是輸出不是壓縮文件汹桦，畢竟壓縮也是需要時間的，因此用內(nèi)存換速度鉴裹，我覺的挺值的

• 修改fastq的名字
  因為CellRanger ATAC只能識別他pipeline規(guī)定的名字格式营勤，因此需要改名，格式如下

[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怎么對應(yīng)I1, R1,R2,R3笆傩摺猖凛？方法很簡單，在瀏覽器中輸入下面的網(wǎng)址查詢一下文件的基本信息绪穆，需要查詢其他的修改run=SRR13450125即可

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

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

• 作者不想讓你分析系列
  如果你探索的足夠深入的話试溯，你還會發(fā)現(xiàn)一些奇怪的東西，比如SRR14539571
  
  
  這種就直接放棄吧郊酒，我看了一下他是把R1(L=51)R2(L=16)R3(L=51)序列合并了遇绞，很難搞，感興趣的同學可以自己把barcode提取出來繼續(xù)分析燎窘，這里不再展開

$ 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系列
  如果還要繼續(xù)探索的摹闽，還存在另外一種情況，就是一個樣本有多個Lane褐健，因此分析的時候要把原始文件下載全了再分析付鹿，就比如說上面下載的SRR13450125，他其實還有個孿生兄弟SRR13450126蚜迅，只不過測序的時候泳道不同舵匾，所以下載的時候一定要下載全，雖然不全也能繼續(xù)分析慢叨，只不過就不能復(fù)原文獻內(nèi)容了纽匙。這種情況在傳統(tǒng)ATAC-Seq中也存在，這里不再展開敘述了

終于把下載數(shù)據(jù)折騰完了

1.3 下載【比對索引】或者【自己構(gòu)建索引】

這里就不得不吐槽10X了拍谐，你整那么大的索引干什么烛缔，下載起來賊費勁馏段。

• 下載適配于Cell Ranger ATAC小鼠的參考基因組（13GB，解壓后19GB）践瓷。沒辦法院喜，就是這么大，但下載總比自己構(gòu)建的快且用的安心吧

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 #快點刪了

• 番外篇：自己構(gòu)建適配于Cell Ranger ATAC的參考基因組——mkref
• 目前官網(wǎng)上提供的參考基因組只有GRCh38與mm10晕翠，因此做其他動物或植物的需要自己構(gòu)建索引喷舀。但必須同時具有g(shù)tf和genome文件，所以不想下載或官網(wǎng)上沒研究的物種的淋肾，也可以自己構(gòu)建

# 首先準備一個配置文件
vi config
# 在配置文件中輸入以下內(nèi)容后保存（字典格式）
{
    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"
}
# 然后執(zhí)行程序
cellranger-atac mkref --config=config

organism：可選參數(shù)硫麻，指定你想定義的物種名
genome：必選參數(shù)，輸出文件所存放的文件夾名稱
input_fasta：必選參數(shù)樊卓，物種的基因組文件
input_gtf：必選參數(shù)拿愧，物種的gtf注釋文件
non_nuclear_contigs：構(gòu)建索引的時候忽略的染色體，chrM指線粒體上的染色體碌尔，指定多個用逗號分割浇辜，如：["chrM","chrC"]
input_motifs:：可選參數(shù)，指定jaspar格式的motif文件唾戚，推薦加上啊柳洋，可以去扣扣群559758504下載hg38或者mm10的motifs文件

1.4 call peak

這里就以我們下載的SRR13450125和SRR13450126為例吧，因為二者是孿生兄弟叹坦，因此分析的時候一起分析熊镣，所以改名的時候不能說改為SRR13450125_S1...和SRR13450126_S1...，要命名成同樣的前綴立由，這里我統(tǒng)一命名為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：下載或自己構(gòu)建的參考基因組索引文件夾
--fastqs：存放樣本的文件夾。路徑前面一定不要有【/】道盏，否則報錯
--sample：樣本名而柑，也就是_S1前面的內(nèi)容，我們前面改名為了armstrong_IGO_09847_1_S1...因此這個地方就是armstrong_IGO_09847_1
--lanes：指定分析一個樣本下的單個lane荷逞，如--lanes=1則只分析L001的內(nèi)容媒咳。如果不使用，則默認分析一個樣本下的所有l(wèi)anes种远，這里我們不使用涩澡，則同時分析L001和L002
--localcores：默認占滿服務(wù)器的所有核CPU，因此推薦設(shè)置的小一點坠敷，如8核
--localmem：請調(diào)整參數(shù)限制占用的運行內(nèi)存

• 注意：想了解更多--fastqs妙同，--sample射富，--lanes如何配合使用，請參考【官網(wǎng)】或【中文解釋】

為了比較差異粥帚，我單獨跑了一下L001和L002的call peak結(jié)果胰耗，差異比較放到結(jié)果解讀后面吧

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 結(jié)果解讀

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!