ixxmu
commented
1 year ago 肝癌单细胞转录组数据之SRP318499 by 生信技能树
这个肝癌单细胞转录组数据的SRA数据库链接是：
https://www.ncbi.nlm.nih.gov/sra/?term=SRP318499 https://trace.ncbi.nlm.nih.gov/Traces/?study=SRP318499
可以看到居然是bam文件格式：
Organism Sample File Name Size Updated
Homo sapiens SAMN19017675 095_possorted_genome_bam.bam 55.2Gb 2021-05-05 10:26:32
Homo sapiens SAMN19017676 104_possorted_genome_bam.bam 22.8Gb 2021-05-05 05:41:14
Homo sapiens SAMN19017677 106_possorted_genome_bam.bam 21.1Gb 2021-05-05 05:51:20
Homo sapiens SAMN19017678 114_possorted_genome_bam.bam 50.1Gb 2021-05-05 08:36:12
Homo sapiens SAMN19017679 119_possorted_genome_bam.bam 44.7Gb 2021-05-05 09:00:40
Homo sapiens SAMN19017680 713_possorted_genome_bam.bam 21.4Gb 2021-05-05 06:10:24
Homo sapiens SAMN19017681 725_possorted_genome_bam.bam 26.2Gb 2021-05-05 05:42:48
Homo sapiens SAMN19017682 740_possorted_genome_bam.bam 52.5Gb 2021-05-05 09:12:26
其实也有对应的bioproject的ID ：https://www.ncbi.nlm.nih.gov/bioproject/PRJNA727404
安装conda
这里推荐每个人安装自己的conda，这样的话一个服务器里面的每个用户独立操作，安装方法代码如下：
# 首先下载文件，20M/S的话需要几秒钟即可
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
# 接下来使用bash命令来运行我们下载的文件，记得是一路yes下去
bash Miniconda3-latest-Linux-x86_64.sh
# 安装成功后需要更新系统环境变量文件
source ~/.bashrc
安装好conda后需要设置镜像（取决于你服务器的情况，部分网络比较好的服务器可以不设置镜像）
conda config --add channels https://mirrors.bfsu.edu.cn/anaconda/cloud/bioconda/
conda config --add channels https://mirrors.bfsu.edu.cn/anaconda/cloud/conda-forge/
conda config --add channels https://mirrors.bfsu.edu.cn/anaconda/pkgs/free/
conda config --add channels https://mirrors.bfsu.edu.cn/anaconda/pkgs/main/
conda config --set show_channel_urls yes
一个完整并且好用的conda真的是生产力工具。接下来就使用conda来配置mamba后设置kingfisher环境。
配置mamba后设置kingfisher环境
七八年前我写教程，尤其是多组学测序数据分析，都是从ncbi的sra数据库里面下载sra文件，然后一步步处理。很多教程还录制成为了视频上传到b站。但是读者多了之后我接受到的大家的反馈就是从ncbi的sra数据库里面下载sra文件实在是太慢了，因为我做演示的服务器在境外，所以自己压根就没有意识到这点。但是陆陆续续有小伙伴告诉我应该是使用aspera从ebi的ena数据库直接下载fastq文件即可，高速而且还少了一个sra文件转为fastq的步骤。所以后来我也开始在日常更新的公众号里面推荐这个方法，就是参考：使用ebi数据库直接下载fastq测序数据 , 需要自行配置好。但是风水轮流转，没想到，也没有过去多少年，就风水轮流转， aspera从ebi的ena数据库这个手段时灵时不灵的。现在只能是回归sra下载，见：转录组上游定量分析其实真不难，4步可定（一），从ncbi的sra数据库里面下载sra文件需要的是sra-tools这个工具套件， 如果你按照我的转录组流程配置：在全新服务器配置转录组测序数据处理环境，会发现安装的是sra-tools-2.8.0 ，后面仍然是报错，所以得指定版本，比如从github安装最新版。
慢慢的，大家在交流群提到了无论是使用sra-tools从ncbi的sra数据库里面下载sra文件，还是使用aspera从ebi的ena数据库直接下载fastq文件，都是得费劲测试和验证， 还得自己准备每个项目的全部的样品的信息，工作繁琐。所以就有了kingfisher，仅仅是需要一个项目ID，就可以自动去检索ebi的ena数据库和ncbi的sra数据库，自动合适的方式下载，非常方便。
代码如下所示：
conda install mamba
mamba create -n kingfisher python=3.8
mamba init
mamba list
mamba activate kingfisher
mamba install -c bioconda kingfisher
一行代码下载数据
mamba activate kingfisher
kingfisher get -r SRP318499 -m ena-ascp ena-ftp prefetch aws-http
因为是根据sra数据库的id所以下载的是sra文件 , 需要把它们转为fq文件
ls *sra| xargs -iF -P 4 sh -c ' kingfisher extract --sra F -t 2 -f fastq.gz'
我看了看，可能是因为作者提供的bam文件的，所以这个过程显得非常奇怪：
ls -lh |cut -d" " -f 5-
13G 1月 20 21:04 SRP318499.1.fastq.gz
13G 1月 20 12:30 SRP318499.1.sra
9.9G 1月 20 21:01 SRP318499.2.fastq.gz
9.5G 1月 20 12:44 SRP318499.2.sra
20G 1月 20 21:49 SRP318499.3.fastq.gz
16G 1月 20 12:59 SRP318499.3.sra
21G 1月 20 21:52 SRP318499.4.fastq.gz
19G 1月 20 13:18 SRP318499.4.sra
19G 1月 20 13:38 SRP318499.5.sra
11G 1月 21 02:21 SRP318499.6.fastq.gz
9.7G 1月 20 13:47 SRP318499.6.sra
12G 1月 21 02:33 SRP318499.7.fastq.gz
9.9G 1月 20 14:02 SRP318499.7.sra
22G 1月 21 03:02 SRP318499.8.fastq.gz
22G 1月 20 14:25 SRP318499.8.sra
289G 1月 20 14:56 SRP318499.fastq
另外，也有可能是因为这个sra数据库项目是10x技术单细胞，所以上面的参数可能是不准确，所以我们重新使用一行代码下载数据，这次根据对应的bioproject的ID ：https://www.ncbi.nlm.nih.gov/bioproject/PRJNA727404
mamba activate kingfisher
kingfisher get -p PRJNA727404 -m ena-ascp ena-ftp prefetch aws-http
但是很不幸，仍然是下载的每个样品是单个fq文件，不符合10x这个单细胞转录组技术。
ls -lh |cut -d" " -f 5-
26G 1月 21 09:54 SRR14424777.fastq.gz
14G 1月 21 10:04 SRR14424778.fastq.gz
11G 1月 21 10:10 SRR14424779.fastq.gz
24G 1月 21 10:25 SRR14424780.fastq.gz
25G 1月 21 10:41 SRR14424781.fastq.gz
12G 1月 21 10:49 SRR14424782.fastq.gz
13G 1月 21 10:57 SRR14424783.fastq.gz
27G 1月 21 11:18 SRR14424784.fastq.gz
这个时候需要参考 只能下载bam文件的10x单细胞转录组项目数据处理 ，如果你熟悉10x单细胞转录组数据，就知道：
首先，1-26个cycle就是测序得到了26个碱基，先是16个Barcode碱基，然后是10个UMI碱基；通常是R1文件
然后，27-34这8个cycle得到了8个碱基，就是i7的sample index；通常是I1文件
最后35-132个cycle得到了98个碱基，就是转录本reads（目前很多测序仪都是150bp了），通常是R2文件
也就是说R2 文件是真正的测序reads，肯定是文件最大。因为这个项目对应的ebi数据库链接是：https://www.ebi.ac.uk/ena/browser/view/PRJNA727404，就是上面的 kingfisher get -p PRJNA727404 命令获取到的8个fq.gz 文件。
所以仍然是得回到 https://trace.ncbi.nlm.nih.gov/Traces/?view=study&acc=SRP318499 页面，去下载那些bam文件，会复杂一点了。需要参考 只能下载bam文件的10x单细胞转录组项目数据处理 进行bam文件转为fq文件，然后再针对fq文件进行cellranger的定量流程。定量后的表达量矩阵很容易进行降维聚类分群，文章给出来的示范如下行：
可以看到，恶性肿瘤上皮细胞是每个病人都不一样，但是其它免疫细胞和间质细胞是可以比较好整合在一起的。感兴趣的可以去读一下文章：
cellranger的定量流程详解：
正常走cellranger的定量流程即可，代码我已经是多次分享了。参考：
10X单细胞转录组原始测序数据的Cell Ranger流程（仅需800元） 10X的单细胞转录组原始数据也可以在EBI下载 一个10x单细胞转录组项目从fastq到细胞亚群 一文打通单细胞上游：从软件部署到上游分析 PRJNA713302这个10x单细胞fastq实战 一次曲折且昂贵的单细胞公共数据获取与上游处理 只能下载bam文件的10x单细胞转录组项目数据处理 不知道10x单细胞转录组样品和fastq文件的对应关系 10X单细胞转录组测序数据的 SRA转fastq踩坑那些事 10x的单细胞转录组fastq文件的R1和R2不能弄混哦
差不多几个小时就可以完成全部的样品的cellranger的定量流程。基础知识非常重要，我们在单细胞天地多次分享过cellranger流程的笔记（2019年5月），大家可以自行前往学习，如下：
单细胞春节福利继续
如果你有自己的单细胞转录组数据需要处理，可以直接参加我们的活动 ，详见：春节期间单细胞转录组数据分析全免费，如果你感兴趣我们处理过的数据集，可以直接填表留下联系方式即可，我们会在元宵节之前统一发放全部的代码和数据给大家的：
文末友情宣传
强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：
生物信息学马拉松授课（买一得五） ，你的生物信息学入门课
https://mp.weixin.qq.com/s/1WwI7kIRVqHZXpqYmfkTnQ