一次错误的单细胞转录组定量

生信技能树

发布于 2024-11-21 09:59:16

12100

代码可运行

运行总次数：0

代码可运行

概述：cellranger是10x的标准上游处理软件，但是只在Linux内核的机器上运行(Ubuntu/Mac)，或者也可以在官网上运行。由于注册原因，尚未成功。 但是官网处理上游文件需要先下载再上传，速度会非常慢，不推荐。首先在SRA网站上下载原始的单细胞测序下机文件

首先是上游文件的下载。在NCBI的GEO上进行搜索，以GSE234052为例。

首先呢，可以看到这个数据集其实是给出来了单细胞表达量矩阵文件，我们仅仅是为了演示上游fq文件处理哈。

可以看到每一个GSM的样本都对应着数个SRA文件或SRR号，这也就意味着每一个样都有好几个测序文件，这可能是因为分了多次测序导致的，每一个样的所有测序文件合并起来才是正确的最终文件。

这里以SRR24809309和SRR24809310的一个样为例进行示范。此处在win系统中只进行查看，下载还是要在Linux中进行。

虚拟机配置 + conda环境配置此处使用了VMware 16和ubuntu 22.04(Jammy Jellyfish 幸运水母version)，centos 7的系统不是很适合，配置起来很麻烦，并且界面不友好。而最新的 ubuntu24.0 对内存的要求似乎比较高，或者是可能会出现掉显卡驱动的事情，导致开机即死机，较难修复。

先下载好ubuntu的iso镜像，注意是桌面版(desktop)不是服务器版，不然在VMware里服务器版装系统时需要比较复杂的配置。新建虚拟机。

在新建虚拟机时，设置好不复杂的用户名和密码，CPU核心数设置推荐2x2最小，多多益善。而硬盘空间的大小可以依需要处理的样本大小而定，选择500GB。内存大小的选用就比较讲究了，在处理共计10GB不到的两个SRA文件时，分配的8GB内存报错，提示程序即需要6GB内存，还需要给其他应用和系统留出内存空间。在分配了20GB内存后解决了这个问题。其他配置基本都不怎么需要修改。具体可以查看VMware安装ubuntu教程。

虚拟机的开机。

成功安装和配置虚拟机之后，就可以开机后可以使用 Ctrl+Alt+t 快捷键打开终端，开始进行conda的环境配置及软件安装。

#安装conda依赖包
> apt install libgl1-mesa-glx libegl1-mesa libxrandr2 libxrandr2 libxss1 libxcursor1 libxcomposite1 libasound2 libxi6 libxtst6

#去linux的浏览器的conda官网找到conda的下载链接，并用 wget 下载
> wget https://repo.anaconda.com/archive/Anaconda3-2022.10-Linux-x86_64.sh

# 安装anaconda
> bash Anaconda3-2020.07-Linux-x86_64.sh

> su #先进入管理员模式给conda配置一次环境变量

# 查看安装脚本配置的环境变量内容
> vim ~/.bashrc

#--文本末尾加入环境变量内容

# >>> conda initialize >>>
# !! Contents within this block are managed by 'conda init' !!
__conda_setup="$('/usr/local/anaconda3/bin/conda' 'shell.bash' 'hook' 2> /dev/null)"
if [ $? -eq 0 ]; then
    eval "$__conda_setup"
else
    if [ -f "/usr/local/anaconda3/etc/profile.d/conda.sh" ]; then
        . "/usr/local/anaconda3/etc/profile.d/conda.sh"
    else
        export PATH="/usr/local/anaconda3/bin:$PATH"
    fi
fi
unset __conda_setup
# <<< conda initialize <<<

#文本编辑器gedit用法
> Esc #退出编辑模式
> i #进入插入模式 进入时默认为插入模式
> :q #无修改退出
> :wq #保存退出
> :q! #不保存修改 强制退出

#输入conda
> exit #退出管理员账户
>conda --help #看看环境变量有没有配置好
#Linux的环境配置非常非常不稳定，可能一直配置不好，也可能重启后又需要重新配置

#创建conda环境
> conda create -n cellranger_09_13

#激活conda环境
> conda activate cellranger_09_13

#添加conda国内外镜像
> conda config --add channels conda-forge   
> conda config --add channels defaults   
> conda config --add channels r   
> conda config --add channels bioconda
> conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
> conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/
> conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/msys2/
> conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/
> conda config --add channels https://mirrors.ustc.edu.cn/anaconda/pkgs/free/
> conda config --add channels https://mirrors.ustc.edu.cn/anaconda/cloud/bioconda/
> conda config --set show_channel_urls yes

#安装sra-tools
> conda install -y -c bioconda sra-tools

之后在Linux浏览器上找到之前的GSE数据集的源文件位置，选择“Accession list”下载，文件中只是简单地写明了需要下载文件的SRR号。

#启用sra-tools的prefetch下载
> prefetch --option-file SRR_Acc_List.txt

#建设好自己的工作文件夹，不要把 文件 和 软件输入输出结果 搞得到处都是
#只提供语法词典，文件夹自定义
> cd ~ #进入home文件夹
> cd - #进入上一文件夹
> cd.. #进入上一层目录
> cd / #进入根目录/计算机
> mkdir workplace #当前目录下创建workplace文件夹
> mv SRR_Acc_List.txt workplace/ #将*.txt文件移到workplace文件夹下

#写一个将*.sra文件批量转为*.fastq格式文件的sh脚本
> touch sra2fastq.sh #创建脚本文件

for i in *sra  
do  
echo $i  
pfastq-dump --gzip --split-files -t 2 $i  #2代表线程数
done

#安装pfastq-dump
git clone https://github.com/inutano/pfastq-dump
cd pfastq-dump
chmod a+x bin/pfastq-dump #给予所有用户使用脚本文件的权限

#将pfastq-dump添加到环境
sudo vim ~/.bashrc
#最后面添加一行,注意修改你保存的绝对路径
export PATH=$PATH:/home/xuran/Desktop/10X/pfastq/pfastq-dump/bin:$PATH
#保存后刷新bash环境
source ~/.bashrc

#给予脚本运行权限并运行
chmod 777 sra2fastq.sh # 777 指全部用户的全部权限
bash sra2fastq.sh #运行

在通过pfastq-dump拆分.sra文件后，可以得到2-4+个*.fastq.gz文件。这些.gz文件都有着不同的作用和意义，可以通过查看测序内容来判定他们的属性。读段比较长的一般应该认为是真正的测序文件，而比较短的片段文件则是索引片段。

> zcat SRR24809309_S1_L001_R1_001.fastq.gz | head -n 10

在确认了各个文件的性质后，应当把测序.gz文件的名字改成可以被cellranger识别的文件名格式。可以看到，[sample name]为了方便识别可以直接命名为SRR号。[S1]一般来讲，同时处理的SRR文件都是来自同一个单细胞样本，同一个单细胞测序样本的各SRR文件应当具有相同的S号。[Lane number]如非特殊关注的情况均默认1即可。[read type]主要分为两种：R1,R2,R3...和I1,I2,I3...可以看到的是，我在之前的工作中已经把文件的名字提前改好了。

这里可以通过sh脚本批量命名，但为了降低复杂度，我们采取手动命名的方式。

> mv old_name.fastq.gz SRR24809309_S1_L001_R1_001.fastq.gz #重命名

之后进行cellranger的处理

#参考基因组下载（human） 大概10个G左右
curl -O https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-GRCh38-2020-A.tar.gz
# mm10
curl -O https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-mm10-2020-A.tar.gz

#cellranger下载，注意要下载最新版本的cellranger 链接会不一样，官网直接找就可以了
curl -o cellranger-8.0.1.tar.gz "https://cf.10xgenomics.com/releases/cell-exp/cellranger-8.0.1.tar.gz?Expires=1726292197&Key-Pair-Id=APKAI7S6A5RYOXBWRPDA&Signature=R9zz3INN5Atmk-69bjlailQwrMPHbZKx3eHLolYmlRKZbZW3e5CK6APP93PQKCtKu6Td2wMNQ-0ay5uGiQyzQtbXKO8gs9JhebHhHe7vE9z3Y645emfnkcruQEPg4v9Izu79JK8BUDOAnrEof77hf5DKcvFAogTHmxZE~aEmr4Bdgj4--sfDA6T~D9W9shFRy2a79F6bCZEjrvXwzQxNEG1n2nKW~dff9Xx7dXYNWKpJ7dJwOx7QVNvjxkfCIPJFGEf6FJWKFLPOX~l3msb7lzsXDdOTKEVifkhMsPozsB3QqxJluySAISRjQtkGZ63RGGYDdMBa8sKa4XFj4Wg-lA__"


#curl -0 #使用http 1.0 
#curl -o #将服务器的回应保存成文件，等同于wget命令
#curl -O #将服务器回应保存成文件，并将 URL 的最后部分当作文件名。

下载好之后解压文件(cellranger + refdata)

> tar -zxvf cellranger-8.0.1.tar.gz
> tar -zxvf refdata-gex-mm10-2020-A.tar.gz

#配置cellranger环境变量
> sudo vim ~/.bashrc #打开环境编辑文件
export PATH=~/download/cellranger-6.1.2:$PATH #加到文件最后一行
> source ~/.bashrc #更新环境变量

> cellranger testrun --id=tiny #测试是否成功
> Pipestance completed successfully! #成功返回值

之后进行cellranger的运行环节，这里为了方便解决报错以及节省内存和精力，我们不选用sh脚本，并手动运行cellranger。

cellranger count --id=SRR24809309 --transcriptome=/study_Linux/workplace/refdata-gex-mm10-2020-A --fastqs=/study_Linux/workplace --sample=SRR24809309 --create-bam=false --nosecondary

#参数注释
cellranger count --id=SRR24809309 \  #id是此次任务名，似乎会影响输出的文件夹名字
--transcriptome=/study_Linux/workplace/refdata-gex-mm10-2020-A \  #参考基因组的位置
--fastqs=/study_Linux/workplace \  #*.fastq.gz文件位置
--sample=SRR24809309 \  #是之前的[sample name]
--create-bam=false \  #是否创建bam格式的文件
--nosecondary  #是否进行分群，因为之后可以在seurat中进行分析，所以选择不分群

#还有其他诸多参数，此处不再展示。运行的不是很慢，大概1-2h运行好
#如果因为修改参数而多次运行cellranger count 可能会出现重复文件的错误，需要删除之前的结果，再重新运行

#还有一组SRR文件
cellranger count --id=SRR24809310 --transcriptome=/study_Linux/workplace/refdata-gex-mm10-2020-A --fastqs=/study_Linux/workplace --sample=SRR24809310 --create-bam=false --nosecondary

全部运行好之后，我们可以在SRR24809309/outs/filtered_feature_bc_matrix中查看输出的三个文件。

这个时候错误产生了

我不知道在哪里看到的教程，说同一个单细胞测序样本得到的两个SRA文件所拆分和分析所得到的文件，现在我们要把他们使用h5文件通过 cellranger aggr 合并起来，如下所示的代码：

> nano aggregation.csv # 创建写明sample和h5文件的csv文件

#列出所有要合并的样本及其路径
sample_id,molecule_h5 #第一行注明属性
SRR24809309,/study_Linux/workplace/SRR24809309/outs/molecule_info.h5 #写好sample和h5文件的位置
SRR24809310,/study_Linux/workplace/SRR24809310/outs/molecule_info.h5

#运行cellranger
cellranger aggr --id=final_files --csv=aggregation.csv
#--id 是输出文件夹的名字

实际上这样拿到的矩阵在进行降维聚类分群的时候是有问题，会得到如下所示的问题UMAP图：