search

ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
120 stars 30 forks source link

一文打通单细胞上游：从软件部署到上游分析 #2169

Closed ixxmu closed 2 years ago

ixxmu commented 2 years ago

https://mp.weixin.qq.com/s/pQyUieRVXohY1enqjN4lLA

github-actions[bot] commented 2 years ago

一文打通单细胞上游：从软件部署到上游分析 by 生信技能树


虽然我分享了：10X单细胞转录组原始测序数据的Cell Ranger流程（仅需800元）以及一个10x单细胞转录组项目从fastq到细胞亚群，但是我发现我一直是从ENA数据库使用aspera进行高速下载，但是它有两个弊端，首先是ena数据库的不稳定，其次是部分数据集的单细胞样品在ena上面并不是R1,R2,I1的3个fastq文件形式。

所以补上这个“生信随笔”的笔记，如果你需要从SRA数据库下载sra文件 ：

下面是公众号“生信随笔”的分享

导语

本文介绍在Linux环境下10X单细胞上游分析环境部署，主要从软件安装和上游分析示例讲起。上游分析示例部分包括两个内容：

  • 常规的SRA文件下载，转fastq格式，然后走cellranger流程；

  • 此外，有些时候作者提供的是bam文件，我们可以先把bam转为fastq，再走cellranger流程。


一. 部署10x上游分析环境

(1) conda创建分析环境及常用软件安装

#conda install mambaconda create -n 10xconda activate 10xmamba install -y -c bioconda aspera-cli bwa samtools bedtools sambamba sra-tools bowtie2 samblaster fasterq-dump

(2) CellRanger及参考基因组下载

2021年12月的最新版本是cellranger软件6.2.1，hg38及mm10参考基因组。

这里cellranger的下载链接会失效，需要去官网获取：https://support.10xgenomics.com/single-cell-gene-expression/software/downloads/latest

注：服务器下载慢的话可以在windows环境下载，然后上传压缩文件至服务器

##1.1 cellranger软件6.2.1wget -O cellranger-6.1.2.tar.gz "https://cf.10xgenomics.com/releases/cell-exp/cellranger-6.1.2.tar.gz?Expires=1636331862&Policy=eyJTdGF0ZW1lbnQiOlt7IlJlc291cmNlIjoiaHR0cHM6Ly9jZi4xMHhnZW5vbWljcy5jb20vcmVsZWFzZXMvY2VsbC1leHAvY2VsbHJhbmdlci02LjEuMi50YXIuZ3oiLCJDb25kaXRpb24iOnsiRGF0ZUxlc3NUaGFuIjp7IkFXUzpFcG9jaFRpbWUiOjE2MzYzMzE4NjJ9fX1dfQ__&Signature=LVCM9NLNJQ4BgDFmIdLG7jYI~hoT4QB4nNIT9U6krn~gVRHQwLP25YDI~msH6cq6fpd8OC4MTNlLNUrqVFTlBK0~qq9Cq~MA-ofUBxYt74yECf4vgGd5uRTUKsiKy2af~6~kicBaj5ltCoIjmguKVQPof4M2E81IdkMUHsVBp6NoCfQCI68sRtIMyYTx7~fpvY1MZv6PipPUXMctB85usluXc~vMQT6f-W4q-gYKG31y8wqS~4Idh4l1oQpMU4T-I16E-EO04~aLQP9lJSMuswZsUxmOP8~wkR93ULmW16UmCH6YYAXAsVZipc6xe1k12yPZ--ioAT9x4JwZW8EIUQ__&Key-Pair-Id=APKAI7S6A5RYOXBWRPDA"
md5sum cellranger-6.1.2.tar#310d4453acacf0eec52e76aded14024c cellranger-6.1.2.tartar -xzvf cellranger-6.1.2.tar
##1.2 mm10wget https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-mm10-2020-A.tar.gz
md5sum refdata-gex-mm10-2020-A.tar.gz#886eeddde8731ffb58552d0bb81f533d refdata-gex-mm10-2020-A.tar.gztar -xzvf refdata-gex-mm10-2020-A.tar.gz
##1.3 hg38wget https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-GRCh38-2020-A.tar.gz
md5sum refdata-gex-GRCh38-2020-A.tar.gz#dfd654de39bff23917471e7fcc7a00cd refdata-gex-GRCh38-2020-A.tar.gztar -xzvf refdata-gex-GRCh38-2020-A.tar.gz
##1.4 环境变量配置：PATH的路径依据自己实际的情况而定ln -s ~/cellrange_soft/cellranger-6.1.2/bin/cellranger ~/anaconda3/envs/10x/bin/
cellranger


二. 常规的单细胞SRR文件格式上游分析

创建分析所用的文件夹备用：

mkdir 1.sra 2.raw_fastq 3.cellranger

(1) 数据下载

示例数据GSE117988：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE117988

Metadata是表型数据，Accession List提供了SRA ID：

cd 1.sra
## 写入需要下载的文件名cat >download_fileSRR7722937SRR7722938SRR7722939SRR7722940SRR7722941SRR7722942
# 批量下载cat download_file |while read id;do (prefetch $id &);done  # 后台下载

(2) fasterq-dump SRA转fastq

常规的SRA转fastq文件，用的是fastq-dump软件，速度非常慢，4-5个小时才能处理完一个样本。

这里用新办法fasterq-dump，2分钟完成一个样本：

cd 2.raw_fastqln -s ../1.sra/SRR* ./
### 方法1：多个文件批量做cat >fastq.shls SRR* | while read id;do ( nohup fasterq-dump -O ./ --split-files -e 6 ./$id  --include-technical & );done## 运行脚本bash fastq.sh
### 方法2：for循环一个一个做for i in `ls SRR*`doi=$iecho "fasterq-dump -O ./ --split-files -e 40 ./$i --include-technical"done >fastq.sh## 运行脚本bash fastq.sh
### 大概耗时2分钟完成ls -lh

如果觉得占空间，可压缩一下：

### 压缩一下ls SRR*fastq | while read id; do gzip $id; done

拿到了这些每个样品的多个fastq文件，后面的流程就跟我分享的：10X单细胞转录组原始测序数据的Cell Ranger流程（仅需800元）以及一个10x单细胞转录组项目从fastq到细胞亚群，是一回事啦。

(3) Cellranger count流程

对10X的fq文件运行CellRangerd的counts流程，先做一个测试：

首先需要对fastq.gz文件改名字，SampleName_S1_L001_R1_001.fastq.gz：

mv SRR7722937_1.fastq.gz SRR7722937_S1_L001_I1_001.fastq.gzmv SRR7722937_2.fastq.gz SRR7722937_S1_L001_R1_001.fastq.gzmv SRR7722937_3.fastq.gz SRR7722937_S1_L001_R2_001.fastq.gZ

再运行cellranger：

cd 3.cellranger
ref=/home/data/vip10t17/software_install/10x_refernce/refdata-gex-GRCh38-2020-Aid=SRR7722937cellranger count --id=$id \--transcriptome=$ref \--fastqs=/home/data/vip10t17/GEO_data/10x_test/2.raw_fastq \--sample=$id \--nosecondary \--localmem=30

会一个的话，就可以写shell脚本批量完成：

##第一步 改fastq文件名字：cd 2.raw_fastq
cat ../1.sar/download_file | 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：ref=/home/data/vip10t17/software_install/10x_refernce/refdata-gex-GRCh38-2020-Als *.fastq.gz | cut -d "_" -f 1 | uniq | while read id;donohup cellranger count --id=$id \--transcriptome=$ref \--fastqs=/home/data/vip10t17/GEO_data/10x_test/2.raw_fastq \--sample=$id \--nosecondary \--localcores=10 \ #设置核心数--localmem=30 &done
##第三步可选 打包压缩文件方便传输tar -zcvf Output.tar.gz SRR7722937 SRR7722938 SRR7722939 SRR7722940 SRR7722941 SRR7722942

最主要的几个参数：

  • --id 指定输出文件夹的名字

  • --transcriptome 指定参考基因组的路径

  • --sample 指定需要处理的fastq文件的前缀

  • --expect-cell 指定预期的细胞数目，默认参数是3000个

  • --localcores 指定计算的核心数

  • --mempercore 指定内存大小 GB

  • --nosecondary 不需要进行降维聚类（因为后期会用R可视化）

  • --chemistry，默认是自动识别chemistry，但是有些时候识别不出chemistry的时候，需要加入参数特别标明

使用cellranger count --help可查看更多参数，Cellranger count结果解读见文末。


三. 如果提供的是BAM数据

注意，bam转fastq，再走cellranger的流程非常耗费计算机资源和时间，八个样本预计会占2T的硬盘空间。练习的话，笔者建议下载1-2个样本即可。

(1) 下载BAM数据：

示例数据来自：https://www.ebi.ac.uk/ena/browser/view/PRJNA727404?show=reads

cd ~/Projects_gao/newhbvhc##3.1构建下载listcat >download_file/vol1/run/SRR144/SRR14424777/740_possorted_genome_bam.bam/vol1/run/SRR144/SRR14424778/725_possorted_genome_bam.bam/vol1/run/SRR144/SRR14424779/713_possorted_genome_bam.bam/vol1/run/SRR144/SRR14424780/119_possorted_genome_bam.bam/vol1/run/SRR144/SRR14424781/114_possorted_genome_bam.bam/vol1/run/SRR144/SRR14424782/106_possorted_genome_bam.bam/vol1/run/SRR144/SRR14424783/104_possorted_genome_bam.bam/vol1/run/SRR144/SRR14424784/095_possorted_genome_bam.bam
##3.2批量下载nohup ascp -v -QT -l 300m -P33001 -k1 -i /home/data/ssy40/anaconda3/envs/10x/etc/asperaweb_id_dsa.openssh --mode recv --host fasp.sra.ebi.ac.uk --user era-fasp --file-list download_file ./ & >download.log
ls -lh | cut -d" " -f 5- 56G 5月   6  2021 095_possorted_genome_bam.bam 23G 5月   9  2021 104_possorted_genome_bam.bam 22G 5月  11  2021 106_possorted_genome_bam.bam 51G 5月  22  2021 114_possorted_genome_bam.bam 45G 5月   6  2021 119_possorted_genome_bam.bam 22G 5月   8  2021 713_possorted_genome_bam.bam 27G 5月   8  2021 725_possorted_genome_bam.bam 53G 12月 11  2021 740_possorted_genome_bam.bam449K 12月 10 13:45 wget-log

(2) bam转fastq

参考官网流程：https://support.10xgenomics.com/docs/bamtofastq?src=pr&lss=none&cnm=&cid=NULL

cellranger产生的bam文件里是带有barcode与UMI的，储存在tag标签里：

samtools view 104_possorted_genome_bam.bam | less -SNsamtools view 104_possorted_genome_bam.bam | head -3 | tr "\t" "\n" | cat -n

  • CB、CR、CY表示barcode，一般是16个碱基；

  • UB、UR、UY表示UMI，一般是10个碱基。

  • R一般代表原始测序数据，Y代表质量分数，而B代表校正后的R，可能对应碱基质量分数太低等因素。一般来说R与B都是相同的。

cellranger bamtofastq:

## 3.3构建一个name list文件cat >name.list740725713119114106104095
##3.4创建文件夹mkdir fastq_file
##3.5准备shell脚本cat > bamtofastq.shcat name.list |while read iddocellranger bamtofastq --nthreads 30 --traceback ${id}_possorted_genome_bam.bam ./fastq_file/${id}done#运行脚本nohup bash bamtofastq.sh & 
#必要时可批量kill任务#ps -ef | grep bamtofastq | awk '{print $2}' | while read id;do kill $id;done  #批量Kill

(3) cellranger count

批量完成：

### 看一下文件夹地址find ~/Projects_gao/newhbvhc/fastq_file/*/*count*/ | grep [1-1000] | grep -v XX/bam 
##3.7 输入剩余未完成的文件listcat >other_file.list/home/data/ssy40/Projects_gao/newhbvhc/fastq_file/095//home/data/ssy40/Projects_gao/newhbvhc/fastq_file/104//home/data/ssy40/Projects_gao/newhbvhc/fastq_file/106//home/data/ssy40/Projects_gao/newhbvhc/fastq_file/114//home/data/ssy40/Projects_gao/newhbvhc/fastq_file/119//home/data/ssy40/Projects_gao/newhbvhc/fastq_file/713//home/data/ssy40/Projects_gao/newhbvhc/fastq_file/725//home/data/ssy40/Projects_gao/newhbvhc/fastq_file/740/
##3.8 批量shell代码cat other_file.list |while read iddoref=/home/data/ssy40/cellrange_soft/10x_refernce/refdata-gex-GRCh38-2020-Asample_name=${id:0-36:3}_out_fileecho "cellranger count --id=$sample_name \--transcriptome=$ref \--fastqs=$id \--sample=bamtofastq \--nosecondary \--localmem=20 \--localcores=30"done>other_file.sh
#运行nohup bash other_file.sh>log_106.114.119.log 2>&1 &


四. Cellranger count结果解读

例如其中一个样本产生的结果：

最重要的结果在outs/文件夹下：

结果解读（参考生信技能树Jimmy的帖子）：

  • web_summary.html：必看，官方说明 summary HTML file ，包括许多QC指标，预估细胞数，比对率等；

  • metrics_summary.csv：CSV格式数据摘要，可以不看；

  • possorted_genome_bam.bam：比对文件，用于可视化比对的reads和重新创建FASTQ文件，可以不看；

  • possorted_genome_bam.bam.bai：索引文件；

  • filtered_gene_bc_matrices：是重要的一个目录，下面又包含了 barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz，是下游Seurat、Scater、Monocle等分析的输入文件，是经过Cell Ranger过滤后构建矩阵所需要的所有文件；

  • filtered_feature_bc_matrix.h5：过滤掉的barcode信息HDF5 format，可以不看；

  • raw_feature_bc_matrix：原始barcode信息，未过滤的可以用于构建矩阵的文件，可以不看；

  • raw_feature_bc_matrix.h5：原始barcode信息HDF5 format，可以不看；

  • analysis：数据分析目录，下面又包含聚类clustering（有graph-based & k-means）、差异分析diffexp、主成分线性降维分析pca、非线性降维tsne，因为我们自己会走Seurat流程，所以不用看；

  • molecule_info.h5：可用于整合多样本，使用cellranger aggr函数；

  • cloupe.cloupe：官方可视化工具Loupe Cell Browser 输入文件，无代码分析的情况下使用，会代码的同学通常用不到。


ixxmu commented 2 years ago

https://share.weiyun.com/R4F8i9Hu

ixxmu commented 2 years ago

5套单细胞数据分析代码，腾讯微云里面：https://share.weiyun.com/R4F8i9Hu 另外，我们团队长期提供单细胞服务，详见：https://mp.weixin.qq.com/s/cvPsPeQOqk52YX3aSoifAQ

如果你不想看这么多，直接学这个单细胞数据集图表复现代码： 链接：https://pan.baidu.com/s/1KmOWq3fkte6JREyyYUxOLw 提取码：tree 有任何疑问都可以发邮件给我，我的邮箱地址是 jmzeng1314@163.com

ixxmu commented 2 years ago

https://mp.weixin.qq.com/s/v8Vidgu1W0C5hoMPAvpFBw