之前都是拿到单细胞counts矩阵进行后续的分析，这次主要基于10X的cellranger软件学习单细胞原始的测序数据（sra/fastq）到表达矩阵的转换。值得注意的是cellranger不会得到最终的表达矩阵，而是如下图所示的每个样本三个配套文件，然后利用Seurat等R包的函数转化为表达矩阵。
笔记首先总结下自己结合具体数据的“顺利”流程，最后再记录下作为新手遇到的坑。
- 数据支持：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE112903（两个样本）
- 对应文献：https://pubmed.ncbi.nlm.nih.gov/29752062/
putout
一、“顺利”的流程
- 插一句，全程使用的linux环境，大家应该都知道~~
1、下载软件
1.1 cellranger
#目前为最新4.0 version
curl -o cellranger-4.0.0.tar.gz "https://cf.10xgenomics.com/releases/cell-exp/cellranger-4.0.0.tar.gz?Expires=1604161511&Policy=eyJTdGF0ZW1lbnQiOlt7IlJlc291cmNlIjoiaHR0cHM6Ly9jZi4xMHhnZW5vbWljcy5jb20vcmVsZWFzZXMvY2VsbC1leHAvY2VsbHJhbmdlci00LjAuMC50YXIuZ3oiLCJDb25kaXRpb24iOnsiRGF0ZUxlc3NUaGFuIjp7IkFXUzpFcG9jaFRpbWUiOjE2MDQxNjE1MTF9fX1dfQ__&Signature=B3B0SLEsIXZRhpMshMZxcZc-3vUERbq8veptEr2RXSsJnffIXjjQzxuZvnQVBFH0TVdP5pZ6vZuirgmYey04tV8AveldFy7hINl~qumuD950z6UXYcHfK1zXHXzn0b2bPnClUVJHzshjcry0hnBLQHFvdjnTETwn0x~KPXJ8II8Dn8je4X~h7u7a1zO65VFHf-iT9M2nAAEwtzvN0cuSiYtpO9yKioH3IWTismVRvHkwTZZPlfSJagZk3pjxq~-VprqDTmOx-KZwqUn~GfG-sE7HcIwpl7UomovfBz0d6DKSiWGhtQCyY0~kcVdUDnAtRioRL8NZvqGTTHY9zgtJXw__&Key-Pair-Id=APKAI7S6A5RYOXBWRPDA"
tar -zxvf cellranger-4.0.0.tar.gz #解压即用
1.2 other biosoft
- 主要可能会用到
aspera-cli与sra-tools两个生信软件，前者用于高速下载原始测序数据；后者主要用于sra到fastq的转换 - 用conda下载生信软件，方便，高效
conda create -n download #by conda & build an environment
conda activate download
conda install -y -c hcc aspera-cli
conda install -y -c bioconda sra-tools
2、下载数据
2.1 下载原始测序数据
- 当然可以使用GEO提供的原始文件链接wget/curl下载，更推荐使用aspera从ebi数据库下载，快很多.
- 首先就要获取测序在EBI数据库的存储链接：https://mp.weixin.qq.com/s/G4UQNUNXqOzeLVypOJIf6Q
- aspera相关教程参考链接
cat > fq.txt #制备链接txt
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR697/008/SRR6976738/SRR6976738.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR697/009/SRR6976739/SRR6976739.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR697/000/SRR6976740/SRR6976740.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR697/001/SRR6976741/SRR6976741.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR697/002/SRR6976742/SRR6976742.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR697/003/SRR6976743/SRR6976743.fastq.gz
cat fq.txt |while read id
do
ascp -QT -l 300m -P33001 \
-i ~/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh \
era-fasp@$id .
done
如上会下载得到6个fastq.gz测序文件，对应关系为
#GSM3090973 uninfected
SRR6976738.fastq.gz 2G 为序列长度为8的index
SRR6976739.fastq.gz 6G 为序列长度为26的barcode(16)+UMI(10)
SRR6976740.fastq.gz 21G 为序列长度为98的转录本reads
#GSM3090974 post-infected
SRR6976741.fastq.gz
SRR6976742.fastq.gz
SRR6976743.fastq.gz
raw fastq
2.2 下载参考基因组数据
- https://support.10xgenomics.com/single-cell-gene-expression/software/downloads/latest?
- 通过GSE链接，了解到实验对象为Mus musculus小鼠
curl -O https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-mm10-2020-A.tar.gz
tar -zxvf refdata-gex-mm10-2020-A.tar.gz
3、修改名字
3.1 修改fastq文件名
- 为了让cellranger能够正确识别出一个样本的三个配套fatsq数据，一般需要规范文件命名。
-
官方给的示例命名方法之一如下
|-- MySample_S1_L001_I1_001.fastq.gz # I1 means the Index
|-- MySample_S1_L001_R1_001.fastq.gz # R1 means the Barcode+UMI
|-- MySample_S1_L001_R2_001.fastq.gz # R2 means the reads
|-- MySample_S1_L002_I1_001.fastq.gz
|-- MySample_S1_L002_R1_001.fastq.gz
|-- MySample_S1_L002_R2_001.fastq.gz
- 按照上面的示例进行手动改名，主要用
mv命令。代码自动改名语句暂时还没有摸清楚
#GSM3090973 uninfected
GSM3090973_S1_L001_I1_001.fastq.gz #SRR6976738
GSM3090973_S1_L001_R1_001.fastq.gz #SRR6976739
GSM3090973_S1_L001_R2_001.fastq.gz #SRR6976740
#GSM3090974 post-infected
GSM3090974_S1_L002_I1_001.fastq.gz #SRR6976741
GSM3090974_S1_L002_R1_001.fastq.gz #SRR6976742
GSM3090974_S1_L002_R2_001.fastq.gz #SRR6976743
3.2 修改fastq header名（optional）
- 如果下次遇到的情况为：每个样本的配套三文件为同一SRR号，就不用执行这一步骤了。
-
由于每个样本的三个fastq文件分别是一个SRR号，即header不同；这样cellranger在执行命令时不能将三者的数据串到一起（别问我是怎么知道的...）
3 SRR
gunzip ../fastq/GSM3090973_S1_L001_I1_001.fastq.gz &
nohup sed 's/6976738/6976740/' GSM3090973_S1_L001_I1_001.fastq > tmp &
nohup gzip tmp &
mv tmp.gz GSM3090973_S1_L001_I1_001.fastq.gz
gunzip ../fastq/GSM3090973_S1_L001_R1_001.fastq.gz &
nohup sed 's/6976739/6976740/' GSM3090973_S1_L001_R1_001.fastq > tmp2 &
nohup gzip tmp2 &
mv tmp2.gz GSM3090973_S1_L001_R1_001.fastq.gz
gunzip ../fastq/GSM3090974_S1_L002_I1_001.fastq.gz &
nohup sed 's/6976741/6976743/' GSM3090973_S1_L002_I1_001.fastq > tmp3 &
nohup gzip tmp3 &
mv tmp3.gz GSM3090974_S1_L002_I1_001.fastq.gz
gunzip ../fastq/GSM3090974_S1_L002_R1_001.fastq.gz &
nohup sed 's/6976742/6976743/' GSM3090973_S1_L002_R1_001.fastq > tmp4 &
nohup gzip tmp4 &
mv tmp4.gz GSM3090974_S1_L002_R1_001.fastq.gz
为什么如此的修改，可参看下面第二部分的踩坑记录
- 综上，就得到了每套三个的fastq的文件名与header一致的结果，可用于后续的分析。
GSM3090973_S1_L001_I1_001.fastq.gz #SRR6976738
GSM3090973_S1_L001_R1_001.fastq.gz #SRR6976739
GSM3090973_S1_L001_R2_001.fastq.gz #SRR6976740
GSM3090974_S1_L002_I1_001.fastq.gz #SRR6976741
GSM3090974_S1_L002_R1_001.fastq.gz #SRR6976742
GSM3090974_S1_L002_R2_001.fastq.gz #SRR6976743
after rename
GSM3090973
4、cellranger
- 这一步就是利用cellranger软件，结合fastq测序文件与参考基因组文件进行比对、count。
- 命令很简单，但参数要求与数据要求都比较严谨的，要注意、小心
cat > GSM3090973.sh
bin=../cellranger-4.0.0/bin/cellranger #设置变量，方便应用
db=./refdata-gex-mm10-2020-A
fq_dir=./fq/
$bin count --id=GSM3090973 \ #输出文件夹名
--localcores=4 \ #线程
--transcriptome=$db \ #参考基因组路径
--fastqs=$fq_dir \ #测序数据路径
--sample=GSM3090973 \ #指定sample的配套文件的前缀，从而在$fq_dir里识别
--expect-cells=3000 \ #可通过文献查询，大致给一个数量
--nosecondary #不进行下游分析
nohup bash GSM3090973.sh 1>GSM3090973.log 2>&1 &
cat > GSM3090974.sh
bin=../cellranger-4.0.0/bin/cellranger
db=./refdata-gex-mm10-2020-A
fq_dir=./fq/
$bin count --id=GSM3090974 \
--localcores=4 \
--transcriptome=$db \
--fastqs=$fq_dir \
--sample=GSM3090974 \
--expect-cells=12000 \
--nosecondary
nohup bash GSM3090974.sh 1>GSM3090974.log 2>&1 &
- 上述两个命令从半夜12点跑到第二天中午12点左右，12h才执行完毕；
- 不过看到最终成功的结果还是挺开心的。
理论上只需进行四部即可达到基本的目标，拿到比对count数据。估计按上述顺利的流程一天一夜就应该能处理好了，可我断断续续花了两天半的时间才完成。作为新手，总会遇到一系列的"坑"
二、踩坑记录
- 因为第一次学习cellranger软件，所以学习过程中碰到了不少坑。大部分通过自行google都可以解决。
- 但也遇到了几个郁闷的问题，导致耗费了不少时间。
1、header mismatch
- 简单来说就是上面3.2步骤解决的问题。
- 一开始未进行3.2的修改，直接运行第四步：提示的报错类似 input data header mismatch之类的报错；google、baidu都没有找到类似的解答，很郁闷，也没想到是fastq的header不一致的问题；
- 于是就想着按照单细胞天地之前的示例教程走一遍，看看到底哪里不一样；按照示例教程最后运行成功了（代码就不演示了，详见原文教程），想了下与上述文件的不同在于：示例中每个样本就提供一个sra文件，利用sam-tools转化为三个fastq，即这三个fastq的名字(SRR ID)是一样的，最多只是后缀不同
- 而本次任务每个样本直接提供三个fastq文件，对应index、barcode+UMI、reads，即有三个SRR号。虽然文件名容易改，但内部的fastq组成的序列header还是不一样的。
2、sra-tools
- 在按照示例教程，将SRA文件转为fastq格式时，一般会用到sra-tools；参照教程，进行如下命令拆分
fastq-dump --split-3 ******sra - 结果提示 --split-3 parameter unrecognized；查看帮助文档也是推荐这个参数设置的，真的是无语了。
- 后来尝试了 --split-files参数成功得拆分成了三个文件
3、change header
- 如第一点所说，发现问题后，就要补坑。其实也不是一个坑，而涉及到效率。同样的目的，如果能花1个小时解决，为什么要花3个小时呢。关于修改这个header换了好几次思路（后面估计cellranger更耗时，就想睡之前把cellranger挂到后台，即尽快把这个change完成）
- try1：将配套的三个fastq都换成统一的name
zless GSM3090973_S1_L001_R2_001.fastq.gz | sed 's/SRR6976738/GSM3090973/' | gzip> ../rename/GSM3090973_S1_L001_R2_001.fastq.gz &
zless GSM3090973_S1_L001_R2_001.fastq.gz | sed 's/SRR6976739/GSM3090973/' | gzip> ../rename/GSM3090973_S1_L001_R2_001.fastq.gz &
zless GSM3090973_S1_L001_R2_001.fastq.gz | sed 's/SRR6976740/GSM3090973/' | gzip> ../rename/GSM3090973_S1_L001_R2_001.fastq.gz &
- try2：后来想第三个fastq文件太大20多G，就把第一、二个fastq文件改成与第三个一致，不是更快？
zless GSM3090973_S1_L001_I1_001.fastq.gz | sed 's/6976738/6976740/' | gzip> ../rename/GSM3090973_S1_L001_I1_001.fastq.gz &
zless GSM3090973_S1_L001_R1_001.fastq.gz | sed 's/6976739/6976740/' | gzip> ../rename/GSM3090973_S1_L001_R1_001.fastq.gz &
#try3
zless ../fastq/GSM3090973_S1_L001_I1_001.fastq.gz | paste - - - - | sed 's/6976738/6976740/' | tr "\t" "\n"| gzip > GSM3090973_S1_L001_I1_001.fastq.gz &
zless ../fastq/GSM3090973_S1_L001_R1_001.fastq.gz | paste - - - - | sed 's/6976739/6976740/' | tr "\t" "\n"| gzip > GSM3090973_S1_L001_R1_001.fastq.gz &
#试了上述两种方法，都还是需要好几个小时
- try4
即第一部分3.2步进行的操作，1-2h即可完成修改。相比上述方法，效率还是很高的。
至于速度差异的原因，目前也不是非常清楚。
