logo
搜索
我的图书馆

发文章

发文工具

其他工具

留言交流
搜索
分享

【原】小鼠的5个样品的10x技术单细胞转录组上游定量

 健明 2023-06-23 发布于广东

新鲜出炉（2023年5月）的文章：《Fueling sentinel node via reshaping cytotoxic T lymphocytes with a flex-patch for post-operative immuno-adjuvant therapy》

看到了一个文献，它的单细胞转录组并没有给出来表达量矩阵，仅仅是测序数据：

The single-cell RNA sequencing data generated in this study have been deposited in the Sequence Read Archive database with the accession code of PRJNA853539.

我就去 https://www.ncbi.nlm./sra?linkname=bioproject_sra_all&from_uid=853539 看了看：

仅仅是测序数据

那就很简单的了，只需要下载这些数据即可，然后走cellranger的定量流程。

cellranger的定量流程详解：

正常走cellranger的定量流程即可，代码我已经是多次分享了。参考：

差不多几个小时就可以完成全部的样品的cellranger的定量流程。基础知识非常重要，我们在单细胞天地多次分享过cellranger流程的笔记（2019年5月），大家可以自行前往学习，如下：

安装conda

自己下载安装自己的conda，自己配置哈。

# 首先下载文件，20M/S的话需要几秒钟即可
wget https://repo./miniconda/Miniconda3-latest-Linux-x86_64.sh
# 接下来使用bash命令来运行我们下载的文件，记得是一路yes下去
bash Miniconda3-latest-Linux-x86_64.sh 
#  安装成功后需要更新系统环境变量文件
source ~/.bashrc

首先如果是在中国大陆，需要设置好镜像：

conda config --add channels https://mirrors./anaconda/cloud/bioconda/
conda config --add channels https://mirrors./anaconda/cloud/conda-forge/
conda config --add channels https://mirrors./anaconda/pkgs/free/
conda config --add channels https://mirrors./anaconda/pkgs/main/
conda config --set show_channel_urls yes 

然后使用自己的conda来安装 kingfisher 软件：

conda install mamba
mamba create -n kingfisher python=3.8
mamba init
mamba list
# 重启terminal：
mamba activate kingfisher
mamba install -c bioconda kingfisher

然后就一行代码下载数据

kingfisher get -p PRJNA853539  -m ena-ascp ena-ftp prefetch aws-http

另外，值得注意的是，因为上面的代码耗时很久，一般来说需要使用nohup包裹一下，或者说开启screen命令：

screen是一个非常有用的命令行工具，它允许你在单个终端窗口中创建、使用和管理多个终端会话。这在你需要在远程服务器上运行长时间的任务或同时运行多个任务时特别有用。

以下是一些基本的screen命令：

  1. screen：启动一个新的screen会话。
  2. screen -S <session_name>：启动一个新的screen会话，并给它一个名字，这样你可以更容易地记住和引用它。
  3. screen -ls：列出当前所有的screen会话。
  4. screen -r <session_name>：恢复（重新连接到）一个已经存在的screen会话。
  5. screen -d <session_name>：断开一个screen会话。这将保持会话中的进程运行，但会话将不再显示在你的终端窗口中。

在一个screen会话中，你可以使用一些特殊的命令来管理你的会话。这些命令通常以Ctrl-a开始，然后跟着一个字符来执行特定的操作。例如：

  1. Ctrl-a c：在当前screen会话中创建一个新的窗口。
  2. Ctrl-a n和Ctrl-a p：在当前screen会话的窗口之间切换。
  3. Ctrl-a d：断开当前screen会话，使其在后台运行。
  4. Ctrl-a k：关闭当前窗口。
  5. Ctrl-a ?：显示screen的帮助信息，列出所有可用的命令和快捷键。

使用screen命令可以让你在单个终端窗口中更有效地管理你的工作，特别是在处理远程服务器或长时间运行的任务时。有类似功能的其他工具：

  1. tmux：这是一个非常强大的终端复用器，它提供了许多screen没有的功能，如窗口分割、会话管理和自定义键绑定。tmux的学习曲线可能比screen稍微陡峭一些，但一旦你习惯了它，你可能会发现它比screen更强大、更灵活。
  2. byobu：这是一个在screen和tmux之上添加了一些额外功能的包装器。它提供了一个友好的用户界面，包括状态栏和系统通知，使得管理你的终端会话更加容易。
  3. dtach：这是一个非常轻量级的工具，它只做一件事：分离和重新连接到终端会话。如果你不需要screen或tmux的所有额外功能，dtach可能是一个好选择。
  4. abduco：这是另一个轻量级的终端复用器，类似于dtach，但提供了一些额外的功能，如会话列表和更好的终端支持。

这些工具都有各自的优点和缺点，所以你可能需要试用一下，看看哪一个最适合你的需求。

查看10x技术单细胞转录组上游测序数据并且定量

前面的kingfisher下载了PRJNA853539里面的小鼠的5个样品的10x技术单细胞转录组上游测序数据，如下所示：

 ls -lh *z|cut -d" " -f 5-
4.5G 6月  21 17:35 SRR19904954_1.fastq.gz
 11G 6月  21 17:47 SRR19904954_2.fastq.gz
4.7G 6月  21 17:50 SRR19904955_1.fastq.gz
 11G 6月  21 17:58 SRR19904955_2.fastq.gz
5.3G 6月  21 18:08 SRR19904956_1.fastq.gz
 13G 6月  21 18:17 SRR19904956_2.fastq.gz
5.2G 6月  21 18:20 SRR19904957_1.fastq.gz
 13G 6月  21 18:26 SRR19904957_2.fastq.gz
4.4G 6月  21 18:28 SRR19904958_1.fastq.gz
 11G 6月  21 18:37 SRR19904958_2.fastq.gz

这些文件名字是需要改名的。。。。如果你熟悉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，肯定是文件最大。而且I1文件是可以省略的。。。

ls *gz|cut -d"_" -f1 |sort -u | while read id ;do   
mv ${id}_1*.gz ${id}_S1_L001_R1_001.fastq.gz;
mv ${id}_2*.gz ${id}_S1_L001_R2_001.fastq.gz;
done

简单的修改名字后，如下所示：

ls -lh *z|cut -d" " -f 5-
4.5G 6月  21 17:35 SRR19904954_S1_L001_R1_001.fastq.gz
 11G 6月  21 17:47 SRR19904954_S1_L001_R2_001.fastq.gz
4.7G 6月  21 17:50 SRR19904955_S1_L001_R1_001.fastq.gz
 11G 6月  21 17:58 SRR19904955_S1_L001_R2_001.fastq.gz
5.3G 6月  21 18:08 SRR19904956_S1_L001_R1_001.fastq.gz
 13G 6月  21 18:17 SRR19904956_S1_L001_R2_001.fastq.gz
5.2G 6月  21 18:20 SRR19904957_S1_L001_R1_001.fastq.gz
 13G 6月  21 18:26 SRR19904957_S1_L001_R2_001.fastq.gz
4.4G 6月  21 18:28 SRR19904958_S1_L001_R1_001.fastq.gz
 11G 6月  21 18:37 SRR19904958_S1_L001_R2_001.fastq.gz

然后配置cellranger的定量流程

正常走cellranger的定量流程即可，代码我已经是多次分享了。参考：

差不多几个小时就可以完成全部的样品的cellranger的定量流程。基础知识非常重要，我们在单细胞天地多次分享过cellranger流程的笔记（2019年5月），大家可以自行前往学习，如下：

理论上应该是去10x的官网下载cellranger软件和数据库文件，不过我们的共享服务器已经给大家准备了，在 /home/data/refdir/10x 文件夹里面

ls -lh  /home/data/refdir/10x |cut -d" " -f 5-

391M 8月   8  2022 cellranger-6.0.1.tar.gz
 11G 8月   7  2022 refdata-gex-GRCh38-2020-A.tar.gz
9.9G 8月   7  2022 refdata-gex-GRCh38-and-mm10-2020-A.tar.gz
9.7G 8月   7  2022 refdata-gex-mm10-2020-A.tar.gz

在Linux中，你可以使用tar命令来解压一个.tar.gz文件，并将解压后的文件放到指定的文件夹。以下是具体的命令：


tar -xzvf file.tar.gz -C /path/to/directory

在这个命令中：

  • x代表解压
  • z代表gzip，表示文件是gzip格式的，需要进行解压
  • v代表verbose，表示显示详细的解压过程
  • f代表file，表示后面跟的是要解压的文件名
  • -C代表change to directory，表示解压后的文件将被放到后面指定的目录

请将file.tar.gz替换为你的文件名，将/path/to/directory替换为你想要放置解压后文件的目录。

那么，我们的命令就是：

mkdir raw run soft 
mv *gz raw
tar -xzvf /home/data/refdir/10x/cellranger-6.0.1.tar.gz  -C ./soft
tar -xzvf /home/data/refdir/10x/refdata-gex-mm10-2020-A.tar.gz  -C ./soft

值得注意的是cellranger是一直在更新的，所以推荐去官网下载最新版，这里cellranger的下载链接会失效 ，需要去官网获取：https://support./single-cell-gene-expression/software/downloads/latest

有了10x的官网下载cellranger软件和数据库文件，以及前面的kingfisher下载了PRJNA853539里面的小鼠的5个样品的10x技术单细胞转录组上游测序数据，接下来就可以跑cellranger啦，因为是共享服务器，如果不着急的话，可以一个个样品慢慢跑代码。也可以5个样品一起运行：

ref=../soft/refdata-gex-mm10-2020-A/
cellranger=../soft/cellranger-7.1.0/cellranger
ls ../raw/*gz|cut -d"_" -f 1 |sort -u|cut -d"/" -f 3 | cut -d "_" -f 1 | uniq | while read id;do
nohup $cellranger count --id=$id \
--transcriptome=$ref \
--fastqs=../raw \
--sample=$id \
--nosecondary \
--localcores=4 \
--localmem=30 &
done

因为是5个样品一起运行，所以每个样品耗时三五个小时，合起来也是四五个小时就足够啦。虽然说5个样品一起运行会很快，但是大家的CPU和内存可能会不够，如果是共享服务器或者云服务器，另外一个门槛是硬盘。。。。

有可能是硬盘也不够哦。。。。如果是这样的话，上面的代码大家千万不要添加 nohup  ，应该是把上面的代码写入脚本后一个个运行。

然后简单的整理一下结果矩阵和报表即可：

pro=outputs
pro_folder=$pro
mkdir -p $pro_folder
ls -lh  $pro_folder
 
mkdir -p $pro_folder/html 
ls */outs/web_summary.html |while read id;do (cp $id $pro_folder/html/${id%%/*}.html );done

mkdir -p $pro_folder/matrix 
ls -d */outs/filtered_feature_bc_matrix |while read id;do (cp -r  $id $pro_folder/matrix/${id%%/*} );done

很容易拿到了小鼠的5个样品的10x技术单细胞转录组上游定量后的表达量矩阵文件和报表，如下所示：

$ tree -h outputs/
[ 128]  outputs/
├── [ 224]  html
│   ├── [2.0M]  SRR19904954.html
│   ├── [2.0M]  SRR19904955.html
│   ├── [2.1M]  SRR19904956.html
│   ├── [2.0M]  SRR19904957.html
│   └── [2.1M]  SRR19904958.html
└── [ 224]  matrix
    ├── [ 160]  SRR19904954
    │   ├── [ 46K]  barcodes.tsv.gz
    │   ├── [284K]  features.tsv.gz
    │   └── [ 73M]  matrix.mtx.gz
    ├── [ 160]  SRR19904955
    │   ├── [ 50K]  barcodes.tsv.gz
    │   ├── [284K]  features.tsv.gz
    │   └── [ 76M]  matrix.mtx.gz
    ├── [ 160]  SRR19904956
    │   ├── [ 47K]  barcodes.tsv.gz
    │   ├── [284K]  features.tsv.gz
    │   └── [ 74M]  matrix.mtx.gz
    ├── [ 160]  SRR19904957
    │   ├── [ 33K]  barcodes.tsv.gz
    │   ├── [284K]  features.tsv.gz
    │   └── [ 59M]  matrix.mtx.gz
    └── [ 160]  SRR19904958
        ├── [ 49K]  barcodes.tsv.gz
        ├── [284K]  features.tsv.gz
        └── [ 82M]  matrix.mtx.gz

全部的表达量矩阵在outputs下面的matrix文件夹，而报表在outputs下面的html文件夹，后续只需要简单的读取即可。大约200行代码需要预先熟悉一下：

## 
### ---------------
###
### Create: Jianming Zeng
### Date:  2023-01-16 20:53:22
### Email: jmzeng1314@163.com
### Blog: http://www./
### Forum:  http://www./thread-1376-1-1.html
### CAFS/SUSTC/Eli Lilly/University of Macau
### Update Log: 2023-01-16  First version 
###
### ---------------


rm(list=ls())
options(stringsAsFactors = F) 
source('scRNA_scripts/lib.R')

###### step1:导入数据 ######   
# 付费环节 800 元人民币
 
dir='outputs/matrix/' 
samples=list.files( dir )
samples

# samples = head(samples,10) 
sceList = lapply(samples,function(pro){ 
  # pro=samples[1] 
  print(pro)  
  sce =CreateSeuratObject(counts =  Read10X(file.path(dir,pro )) ,
                          project = gsub('^GSM[0-9]*_','',pro)  ,# pro, #
                          min.cells = 5,
                          min.features = 500 )
  return(sce)
})
names(sceList) 

# gsub('^GSM[0-9]*','',samples)
sce.all=merge(x=sceList[[1]],
              y=sceList[ -1 ],
              add.cell.ids =  gsub('^GSM[0-9]*_','',samples)    )
# gsub('_gene_cell_exprs_table.txt.gz','',
#      gsub('^GSM[0-9]*_','',samples) )  

as.data.frame(sce.all@assays$RNA@counts[1:10, 1:2])
head(sce.all@meta.data, 10)
table(sce.all$orig.ident) 


# 分组信息，通常是隐含在文件名，样品名字里面
# phe=str_split( colnames(sce.all),'[-_]',simplify = T)
# head(phe)
# table(phe[,1])
# sce.all$group= phe[,1]

#sce.all$group=toupper( substring(colnames(sce.all),1,5))
# sce.all@meta.data$group =  ifelse(grepl('ne',sce.all$orig.ident),'ne','e')
# sce.all$group=  gsub( '^[0-9]*_' ,'',sce.all$orig.ident) 
# sce.all$group = sce.all$orig.ident
#table(sce.all@meta.data$group)
table(sce.all@meta.data$orig.ident)


###### step2:QC质控 ######
dir.create("./1-QC")
setwd("./1-QC")
# 如果过滤的太狠，就需要去修改这个过滤代码
source('../scRNA_scripts/qc.R')
sce.all.filt = basic_qc(sce.all)
print(dim(sce.all))
print(dim(sce.all.filt))
setwd('../')


sp='mouse'

###### step3: harmony整合多个单细胞样品 ######
dir.create("2-harmony")
getwd()
setwd("2-harmony")
source('../scRNA_scripts/harmony.R')
# 默认 ScaleData 没有添加"nCount_RNA", "nFeature_RNA"
# 默认的代码都是可以自由修改的，需要对全部的代码有基本的认识。
sce.all.int = run_harmony(sce.all.filt)
setwd('../')

###### step4: 降维聚类分群和看标记基因库 ######
# 原则上分辨率是需要自己肉眼判断，取决于个人经验
# 为了省力，我们直接看  0.1和0.8即可
table(Idents(sce.all.int))
table(sce.all.int$seurat_clusters)
table(sce.all.int$RNA_snn_res.0.1) 
table(sce.all.int$RNA_snn_res.0.8) 

getwd()
dir.create('check-by-0.1')
setwd('check-by-0.1')
sel.clust = "RNA_snn_res.0.1"
sce.all.int <- SetIdent(sce.all.int, value = sel.clust)
table(sce.all.int@active.ident) 

source('../scRNA_scripts/check-all-markers.R')
setwd('../') 
getwd()

dir.create('check-by-0.8')
setwd('check-by-0.8')
sel.clust = "RNA_snn_res.0.8"
sce.all.int <- SetIdent(sce.all.int, value = sel.clust)
table(sce.all.int@active.ident) 
source('../scRNA_scripts/check-all-markers.R')
setwd('../') 
getwd()

last_markers_to_check


###### step5: 确定单细胞亚群生物学名字 ######
# 一般来说，为了节省工作量，我们选择0.1的分辨率进行命名
# 因为命名这个步骤是纯人工 操作
# 除非0.1确实分群太粗狂了，我们就选择0.8 
source('scRNA_scripts/lib.R')
# sce.all.int = readRDS('2-harmony/sce.all_int.rds')
# sp='mouse'
colnames(sce.all.int@meta.data) 
table(sce.all.int$RNA_snn_res.0.8)

# keratinocyte：SFN
#melanocyte1 "MLANA" 
selected_genes=c('MUC4', 'PI3', 'SIX3', # nose
                 'SCGB1A1', 'TFF3',   'KRT1',   'KRT10',    
                 'KRT5', 'TP63',   'DLK2',
                 'MKI67', 'TOP2A', 'CDC20' ,
                 'DMD','PUS7','PGAP1','IFI44L',"MLANA" ,'SFN',
                 'MUC5AC' , 'MUC5B',# secretory cells
                 'FOXJ1', 'TPPP3',  'SNTN',  # multiciliated cells 
                 'DEUP1', 'FOXN4',  'CDC20B' # deuterosomal cells
)
# p1 <- DotPlot(sce.all.int, features = unique(str_to_upper(selected_genes)),
#               assay='RNA' ,group.by = 'RNA_snn_res.0.8'  )  + coord_flip() # +th
# 
# p1

gplots::balloonplot(table(sce.all.int$RNA_snn_res.0.8,sce.all.int$orig.ident))

# 前面的出图需要仔细看
# 理论上可以把细胞周期回归掉
# 也可以把文库大小回归掉

# 付费环节 800 元人民币
# 纯手工操作，非常耗时耗力
# 这个时候可以自由选择不同的分辨率
# 一般来说，为了偷懒就选择0.1的分辨率，但是不建议
# 如果是自己的项目，建议起码选择0.8分辨率以上。
# 会有不同的结果，但是影响不大。

if(F){
  sce.all.int
  celltype=data.frame(ClusterID=0:10 ,
                      celltype= 0:10) 
  #定义细胞亚群        
  celltype[celltype$ClusterID %in% c( 0,5 ),2]='Bcells' 
  celltype[celltype$ClusterID %in% c( 4 ),2]='plasma'   
  celltype[celltype$ClusterID %in% c( 6,8,10 ),2]='myeloids'   
  celltype[celltype$ClusterID %in% c( 1,9 ),2]='cd4_conv_Tcells' 
  celltype[celltype$ClusterID %in% c( 2 ),2]='cd8_Tcells'   
  celltype[celltype$ClusterID %in% c( 3 ),2]='cd4_Treg_Tcells' 
  celltype[celltype$ClusterID %in% c( 7 ),2]='NK'     
  
  head(celltype)
  celltype
  table(celltype$celltype)
  sce.all.int@meta.data$celltype = "NA"
  
  for(i in 1:nrow(celltype)){
    sce.all.int@meta.data[which(sce.all.int@meta.data$RNA_snn_res.0.1 == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]}
  Idents(sce.all.int)=sce.all.int$celltype
  
  sel.clust = "celltype"
  sce.all.int <- SetIdent(sce.all.int, value = sel.clust)
  table(sce.all.int@active.ident) 
  
  dir.create('check-by-celltype')
  setwd('check-by-celltype')
  source('../scRNA_scripts/check-all-markers.R')
  setwd('../') 
  getwd()
}

phe=sce.all.int@meta.data
save(phe,file = 'phe.Rdata') 

这个代码并不是最终的解决方案，仅仅是方便快速查看一个单细胞转录组项目而已，因为就耗时十分钟不到就可以出几十张图表，很方便的了解一个单细胞转录组项目的基本情况。而且我们也使用了这个代码处理了成百上千的公共数据集。可以看到基本上是就跟文章一模一样了:

基本上是就跟文章一模一样了

有了这样的快速认知，后续就可以很容易的去深度解析这个数据集，进而复用它。如果你也想完成这个文章的单细胞环节的复现：《Fueling sentinel node via reshaping cytotoxic T lymphocytes with a flex-patch for post-operative immuno-adjuvant therapy》小鼠的5个样品的10x技术单细胞转录组上游定量，你首先需要服务器，如果你比较拮据，大概率是没办法自己采购自己的实体机服务器

    转藏 分享 献花（0）

    来自： 健明 > 《待分类》

    举报/认领

    0条评论

    发表

    请遵守用户 评论公约

    类似文章