单细胞学习2：10xgenomics从头开始分析

2020-08-10 本文已影响0人 wo_monic

[参考链接](https://www.jianshu.com/p/705a348e8623
单细胞实战1,2,3,4,5 https://my.oschina.net/u/4503882/blog/4438601

1.下载数据
使用prefetch或者ascp都可以。
EBI数据库地址
获取下载命令

ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/srr/SRR772/002/SRR7722942 ./
ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/srr/SRR772/001/SRR7722941 ./
ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/srr/SRR772/000/SRR7722940 ./
ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/srr/SRR772/009/SRR7722939 ./
ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/srr/SRR772/008/SRR7722938 ./
ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/srr/SRR772/007/SRR7722937 ./

我们服务器有时候，会报错,查询解决方案说是33001端口没开，但是过一会儿再试，又能正常下载。

cat >SRR_Acc_List-2586-4.txt
SRR7722937
SRR7722938
SRR7722939
SRR7722940
SRR7722941
SRR7722942

2.安装软件cellranger 4.0版（截止2020.8.11）
下载地址https://support.10xgenomics.com/single-cell-gene-expression/software/downloads/latest 自己注册后，会给出下载地址
下载cellranger软件

wget -O cellranger-4.0.0.tar.gz "https://cf.10xgenomics.com/releases/cell-exp/cellranger-4.0.0.tar.gz?Expires=1597153912&Policy=eyJTdGF0ZW1lbnQiOlt7IlJlc291cmNlIjoiaHR0cHM6Ly9jZi4xMHhnZW5vbWljcy5jb20vcmVsZWFzZXMvY2VsbC1leHAvY2VsbHJhbmdlci00LjAuMC50YXIuZ3oiLCJDb25kaXRpb24iOnsiRGF0ZUxlc3NUaGFuIjp7IkFXUzpFcG9jaFRpbWUiOjE1OTcxNTM5MTJ9fX1dfQ__&Signature=Xx0erSX~dD-bjM4iSNM-PbDrdP2Q~QPNw0qEnyB~hff7U243dKjoeWyukXoHpddkqYDdQ0-HHMOITdPA4CmNXVbGXlMGsnY-LSBySee7mgy~MhXEGBiwDOhTLhoxOgyviZfyW-eMYuutVW6Hy-I3tQjqD1bAENfbVNekD25RkNgkTEsOwkhZzLU1hGRaRYD0kMD10~KCUN7-u0OO9WdKsEw7tzPCvb7Hz85O08ENNabAvnyPLw4ZskHNMmP-ciFTi7IItTqLFIQw8D~HAN5GDsTv14Rd2LfOQnS6OwLWffu-eew9m1FzH7zRzjWgb1AFFwDuXx5MHI~knkSYGYCpsQ__&Key-Pair-Id=APKAI7S6A5RYOXBWRPDA"

下载配套的人类基因组11G

wget -c https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-GRCh38-2020-A.tar.gz

下载小鼠基因组

wget -c https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-mm10-2020-A.tar.gz

cellanger安装方法：

tar -zxvf cellranger-4.0.0.tar.gz
cd cellranger-4.0.0
pwd #获取目录
echo `PATH=上述目录\bin:$PATH` >>~/.bash_profile
source export
cellranger #测试命令即可运行

数据处理
3.1 SRR数据解压

##1解压缩文件
cat SRR_Acc_List-2586-4.txt|while read id
do
    fastq-dump --gzip --split-files ./0.raw/$id
done

#2文件批量重命名
# 将原来的_1.fastq.gz改成_S1_L001_I1_001.fastq.gz
# 依次类推，将原来_2的改成R1，将_3改成R2。批量处理，就利用下载SRA的SRR ID好了。
cat SRR_Acc_List-2586-4.txt | while read i 
do
 mv ${i}_1*.gz ${i}_S1_L001_I1_001.fastq.gz;
 mv ${i}_2*.gz ${i}_S1_L001_R1_001.fastq.gz;
 mv ${i}_3*.gz ${i}_S1_L001_R2_001.fastq.gz;
done

一定要按照上面的要求重命名，不然后面cellranger分析的时候，不会自动识别文件名。

选做分析QC

#质控QC
###config里面存的是重命名后的fastq.gz文件的文件名。
cat config | while read id
do
        fastqc --outdir ./2.qc/ --threads 24 ./1.fastq/${id} >> ./2.qc/${id}_fastqc.log 2>&1
done
multiqc  ./2.qc/*zip  -o ./2.qc/multiqc

multiqc目录有最终的html文件。可以看到各组样本的数据质量。
这批数据，个别有问题。

3.cellranger分为4个pipeline
3.1 cellranger mkfastq
将原始碱基（BCL）文件整合为fastq格式，实质是封装的Illumina bcl2fastq
3.2 cellranger count
将获取的fastq文件比对、过滤、计数，这是一个自动化过程
3.3 cellranger aggr
将多个测序结果归一化到相同的测序深度，重新计算feature-barcode matrices
3.4 cellranger reanalyze

先运行测试文件，测试cellranger能否正常运行。

cellranger testrun --id=testruncellranger id指定运行后保存的目录。
分配24个cpu,需要647s.差不多10分钟多。
查看标准输出，标准输出文件最后出现如下，代表测试通过。

Pipestance completed successfully!

2020-08-11 17:19:13 Shutting down.

下载好参考基因组后，解压缩。

tar -zxvf refdata-gex-GRCh38-2020-A.tar.gz
如果想自己构建参考基因组，是整个流程的难点。
3.2 使用cellranger count开始分析

####参数说明：
#id指定输出文件存放目录名称
#转录组指定与CellRanger兼容的参考基因组
#fastqs指定mkfastq或自定义的fastq文件放置位置
#sample要和fastq文件的初始中的样本保持一致，作为软件识别的标志
#except-cells指定复现的细胞数量，这个要和实验设计结合起来
#nosecondary仅获得表达矩阵，不进行后续的降维，聚类和可视化分析（因为后续会自行用R包去做）

##输如参数时，一定不要有多余的空格，不然会报错。
cat SRR_Acc_List-2586-4.txt|while read id
do
 cellranger count --id=$id --transcriptome=~/soft/cellranger/refdata-gex-GRCh38-2020-A --fastqs=1.fastq --sample=$id --expect-cells=5000 --nosecondary
done

这里没有指定使用的内存和cpu数量，大服务器，内存和cpu足够大，所以就没指定，如果内存小于64G,请指定内存。

3.3 合并上面的分析

3.3.1构建Aggregation CSV

AGG123_libraries.csv

library_id,molecule_h5
LV123,/opt/runs/LV123/outs/molecule_info.h5
LB456,/opt/runs/LB456/outs/molecule_info.h5
LP789,/opt/runs/LP789/outs/molecule_info.h5

其中molecule_h5：文件molecule_info.h5 file的路径

echo library_id,molecule_h5 >AGG123_libraries.csv
cat SRR_Acc_List-2586-4.txt|while read id
do
echo $id,$id/outs/molecule_info.h5 >>AGG123_libraries.csv
done

3.3.2 合并

cellranger aggr --id=AGG123 --csv=AGG123_libraries.csv --normalize=mapped

结果输出到AGG123这个目录中

在AGG123的输出结果中，

newplot (1).png

newplot.png

其中有1个红色的alert，Fraction of Reads Kept (SRR7722940) 16.0%的深度不足50%。

4. 使用Seurat分析（其实cellranger也可以自动分析）

5.难点解决，基因组的构建

论文原文是使用人类基因组和MCPyV病毒基因组一起分析的。