单细胞上游分析

从原始数据——得到文件夹filtered_gene_bc_matrices【其中包含barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz，是下游Seurat、Scater、Monocle等分析的输入文件。】

原始数据下载

现在大多数测序公司以及公共数据库给的都是上面所指的三个表达矩阵，不需要下载原始数据处理。但是如果我们希望用公共数据做一些特殊分析（比如lncRNA,repeat elements）就要重头处理数据了，因为用到的注释文件可能不同

# 安装软件
wget https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.11.1/sratoolkit.2.11.1-centos_linux64.tar.gz
tar zvxf sratoolkit.2.11.1-centos_linux64.tar.gz
# 添加到环境变量
# 添加环境变量
echo 'export PATH=$PATH:/LY_software/sratoolkit.2.11.1-centos_linux64/bin' >> ~/.bashrc
source  ~/.bashrc

到GEO数据库中找GSE117988，获取accession list，这里有两种方法推荐第二种。

SRP/SRA Run Selector

  1. 点击SRP，send to file可以下载
  2. 点击SRA Run Selector可以直接获得Accession List，并且这个界面可以查看Metadata等更详细的信息。点击每个run进入Run Browser——Data access，这里一定要仔细看，有的原始文件是fastq，但有些是bam，bam文件后续需要用bamtofastq。这个软件解压后是二进制文件（binary），一般这种文件都存在初始权限问题，修改权限后方可运行。
# 下载软件
wget https://github.com/10XGenomics/bamtofastq/releases/download/v1.3.5/bamtofastq_linux
# 修改权限
chmod 700 bamtofastq_linux
# 写入路径
echo 'export PATH=/home/biosoftware/bamtofastq_linux:$PATH' >>~/.bashrc
source ~/.bashrc
# 运行时必须要输入该文件名的全程
[LY@logina1 cellranger-tiny-bcl-1.2.0]$ bamtofastq_linux
bamtofastq v1.3.5
Invalid arguments.

Usage:
  bamtofastq [options] <bam> <output-path>
  bamtofastq (-h | --help)


# 将下载的accession_list上传到服务器，或者复制粘贴写入新文件
cat > SRR_Acc_List-2586-4.txt <<END
>SRR7722937
>SRR7722939
>SRR7722941
>SRR7722938
>SRR7722940
>SRR7722942
>END
# 查看文件
cat SRR_Acc_List-2586-4.txt
SRR7722937
SRR7722939
SRR7722941
SRR7722938
SRR7722940
SRR7722942
# &&的含义https://blog.csdn.net/chinabestchina/article/details/72686002
## -O|--output-directory <directory>所以'pwd'指的并不是当前目录，而是文件夹的名字
### 这里pwd带引号可能是怕被当成命令
 [ly@logina1 raw]$ cat SRR_Acc_List-2586-4.txt | while read i
> do prefetch $i -O 'pwd' && echo  "**${i}.sra done**"
> done
# cat SRR_Acc_list.txt |while read i; do prefetch $i -O 'raw' && echo *${i}.sra done**"

生成fastq文件

10x测序结果通常有3个fastq文件：
PBMC_DiscAR_S5_L002_I1_001.fastq.gz
PBMC_DiscAR_S5_L002_R1_001.fastq.gz
PBMC_DiscAR_S5_L002_R2_001.fastq.gz

测序最初的Base calling结果是bcl文件，cellranger软件的mkfastq 命令通过调用 Illumina’s bcl2fastq生成3个标准的fastq文件
官网教程：Cell Ranger mkfastq -Software -Single Cell Gene Expression -Official 10x Genomics Support

# 查看官网cellranger mkfastq指导给出的bcl文件
[Ly@logina1 cellranger-tiny-bcl-1.2.0]$ tree Data/
Data/
└── Intensities
    ├── BaseCalls
    │   └── L001
    │       ├── C100.1
    │       │   └── s_1_1101.bcl.gz
    │       ├── C10.1
    │       │   └── s_1_1101.bcl.gz
    │       ├── C101.1
    │       │   └── s_1_1101.bcl.gz
    │       ├── C102.1
    │       │   └── s_1_1101.bcl.gz

如果环境中没有bcl2fastq则无法运行bcl2fastq，需要安装

# 在官网上注册并下载
wget https://files.softwaredownloads.illumina.com/e8ed3335-5201-48ff-a2bc-db4bfb792c85/bcl2fastq2-v2-20-0-linux-x86-64.zip?Expires=1630143759&Policy=eyJTdGF0ZW1lbnQiOlt7IlJlc291cmNlIjoiaHR0cHM6Ly9maWxlcy5zb2Z0d2FyZWRvd25sb2Fkcy5pbGx1bWluYS5jb20vZThlZDMzMzUtNTIwMS00OGZmLWEyYmMtZGI0YmZiNzkyYzg1L2JjbDJmYXN0cTItdjItMjAtMC1saW51eC14ODYtNjQuemlwIiwiQ29uZGl0aW9uIjp7IkRhdGVMZXNzVGhhbiI6eyJBV1M6RXBvY2hUaW1lIjoxNjMwMTQzNzU5fX19XX0_&Signature=C3MKu~2Of9hWKcl1C-rD8tKv6Iz07lSlSp4mxb-p4cAzOqZi3mKJX9K02T2HvwlmRzt0cZkuB2AQqgKAK~QemZSOqFsHcdAa~EDszuooXqdPd~JuS8qGcJxiorI4hQ26sdJoLoHl1MbZ0gvjCIRZOkFrJRPjmDiYsF1To2-u~ZTluL8ppMwUzUvfy9VRpQHgySKaYmZ24md2xAT~1ffrjvd4ile6rrZspMaKFST4BwuBwItGp3XTpStoQhyaVXYIexETmzYf9OS84K-YwHWBJ5ARSviF0UTfmwW2pgwS7qdOwLVb9pbELVLhojQq~J1zFYXzQGJxH39r5TGaEguAww__&Key-Pair-Id=APKAJO3UYWPXK4A26FPQ
# 解压得到.rpm文件
unzip bcl2fastq2-v2-20-0-linux-x86-64.zip
# 破解.rpm文件
rpm2cpio ./bcl2fastq2-v2.20.0.422-Linux-x86_64.rpm | cpio -idmv
# 得到文件夹
./usr/local/bin/bcl2fastq
# 修改权限chmod 700···
# 单独运行软件会报错，无需理会，cellranger可以正常调用
[Ly@logina1 cellranger-tiny-bcl-1.2.0]$ bcl2fastq
BCL to FASTQ file converter
bcl2fastq v2.20.0.422
Copyright (c) 2007-2017 Illumina, Inc.
···
Error>std::exception::what: Unable to find filter file for lane: 1 and tile: 1101
# 写入路径
echo 'export PATH=/LY_software/usr/local/bin:$PATH' >> ~/.bashrc
source  ~/.bashrc

# 查看cellranger mkfastq 参数
[Ly@logina1 LY_software]$cellranger mkfastq --help
Required:
    --run=PATH          Path of Illumina BCL run folder.

 --id=NAME           Name of the folder created by mkfastq. If not supplied,will default to the name of the flowcell referred toby the --run argument.
 --csv=PATH
# --localmem Restricts cellranger to use specified amount of memory (in GB) to execute pipeline stages. By default, cellranger will use 90% of the memory available on your system. 避免cellranger耗尽所有资源 

# 查看一下文件的内容
[Ly@logina1 LY_software]$ cat exercise/Data/cellranger-tiny-bcl-simple-1.2.0.csv
Lane,Sample,Index
1,test_sample,SI-P03-C9
# 运行
cellranger mkfastq --id=tiny_test --run=Data/cellranger-tiny-bcl-1.2.0 --csv=Data/cellranger-tiny-bcl-simple-1.2.0.csv --localmem=10

查看得到的结果

# 我们想要的三个文件在tiny_test/outs/fastq_path/H35KCBCXY、test_sample下
[Ly@logina1 exercise]$ tree tiny_test/outs/fastq_path/
tiny_test/outs/fastq_path/
├── H35KCBCXY
│   └── test_sample
│       ├── test_sample_S1_L001_I1_001.fastq.gz
│       ├── test_sample_S1_L001_R1_001.fastq.gz
│       └── test_sample_S1_L001_R2_001.fastq.gz
├── Reports
···

结果的命名形式与之前提供的csv文件有关


cellranger-tiny-bcl-simple-1.2.0.csv

需要理解一下这三个文件的构成：
I1：index
R1：read 1
R2：read 2

[Ly@loginb2 test_sample]$ zcat test_sample_S1_L001_I1_001.fastq.gz|less -S

@D00547:905:H35KCBCXY:1:1101:19188:87078 1:N:0:AGATCGGG
AGATCGGG
+
.<<....<
@D00547:905:H35KCBCXY:1:1101:14415:42805 1:N:0:CGTCGTCG
CGTCGTCG
+
..<.....
@D00547:905:H35KCBCXY:1:1101:17646:100743 1:N:0:CATCGTTG
CATCGTTG
+
.......<
@D00547:905:H35KCBCXY:1:1101:13481:41873 1:N:0:GTCATACA
GTCATACA
...

I1的每行都是8nt

[Ly@loginb2 test_sample]$ zcat test_sample_S1_L001_R1_001.fastq.gz|less -S
@D00547:905:H35KCBCXY:1:1101:19188:87078 1:N:0:AGATCGGG
CGCTTCGGCGTTATAACCTCACACTC
+
GGGGGIIIIIIIIIIGGGGIGGGIGG
@D00547:905:H35KCBCXY:1:1101:14415:42805 1:N:0:CGTCGTCG
TCAAACGGCCTGTCTCATCATGGAAG
+
GGGGGIIIIIIIIIIIIIIIIIIIII
@D00547:905:H35KCBCXY:1:1101:17646:100743 1:N:0:CATCGTTG
ACCAATGACCAAATCAAAGAAATGAC
+
GGGGGGIIIIIGGIIIIIIIIIIIII
@D00547:905:H35KCBCXY:1:1101:13481:41873 1:N:0:GTCATACA
AACGTTGGTTAAAGTGTCACAAAGGT
+
GGGGGIIIIIIIIIIIIIIIIIIIII
@D00547:905:H35KCBCXY:1:1101:17238:58729 1:N:0:GTCCTATA
ACTAAAAACCAAGCTGTCGCTACTTC
+
...

R1每行都是26nt

[Ly@loginb2 test_sample]$ zcat test_sample_S1_L001_R2_001.fastq.gz|less -S

@D00547:905:H35KCBCXY:1:1101:19188:87078 2:N:0:AGATCGGG
TCCGTTCTGGTGATTCGTCTAAGAAGTTTAAGATTGCTGAGGGTCAGTGGTATCGTTATGCGCCTTCGTATGTTTCTCCTGCTT
+
GAGGGGGIIGIG<<GGGIGGGGGAGGGGGIGIGGIIIIIGGIGGGGGGGGIGIIGGGIIGGGGIIIGIIIGIIGIIIIGIIIGG
@D00547:905:H35KCBCXY:1:1101:14415:42805 2:N:0:CGTCGTCG
CATGCGGCATACGCTCGGCGCCAGTTTGAATATTAGACATAATTTATCCTCAAGTAAGGGGCCGAAGCCCCTGCAATTAAAATT
+
GGGAAIIGAGGGGIGIGIIIIGIIIGIGGIGIIGIIGIIGGGIIIIIIIIIIIGIIIIGGGGGIGIIIIIIIIIIIGGGGIIIG
@D00547:905:H35KCBCXY:1:1101:17646:100743 2:N:0:CATCGTTG
TAAGAGCCTCGATACGCTCAAAGTCAAAATAATCAGCGTGACATTCAGAAGGGTAATAAGAACGAACCATAAAAAAGCCTCCAA
+
GGGAGAGGGIAGGGGGGGGGIGGGIIGIIIIGIIIIIIIGIGGGGIIIGIGGGGGGAGIIIIIIIGIIIGIGIGIGIIGGIGIA
@D00547:905:H35KCBCXY:1:1101:13481:41873 2:N:0:GTCATACA
CTCACAGTGCGCCGCCATGTGCCTGCAGTGTGGGTGCTGCTCAGCCGGGACCCCCTGGACCCCAATGAGTGTGGTTACCAACCC
+
AGGAGGIGGGAGGAGAGGGGGGGIIGGGIGGGGGGGGIGIIGGGGGGGA<GGIIIIIIIIIIGIIIIIIIIIIIIIGAGGGGGG
@D00547:905:H35KCBCXY:1:1101:17238:58729 2:N:0:GTCCTATA
AGAAATATCCGAAAGTGTTAACTTCTGCGTCATGGAAGCGATAAAACTCTGCAGGTTGGATACGCCAATCATTTTTATCGAAGC
+
GAGAAAGGGGG.<GAGGGGGGGGGGIGGGGIIIIIIIIIIGIIIIIIIIIIIIIIIIIGGGGIGIGIIIIIIIIGIIIIGGGGI
@D00547:905:H35KCBCXY:1:1101:15513:65705 2:N:0:GTCATATG
CTTATGCTAATTTGCATACTGACCAAGAACGTGATTACTTCATGCAGCGTTACCATGATGTTATTTCTTCATTTGGAGGTAAAA
+
GAAAAGIGGAGGGGIGGGGIIIIGGIIGGGIIIIIIIIIIIIIIIIIIGIIIIGIIIIIIIIIIGIGIIIIIIIGGGGIIIGGI
...

R2每行大概84nt

利用Sratool的fastq-dump 同样可以得到三个文件，但是需要重新命名。
认识I1/R1/R2的内容需要详细了解10x的测序原理以及引物设计

R1R2

read1 1:28 / read2 2:98 代表R1和R2 的长度区间
10x Genomics RNA-seq 原理详解 - 简书 (jianshu.com)
R1和R2是什么
cell ranger分析流程

cellranger count得到下游分析文件

这里有测试Demo，官网提供人和鼠的参考基因组以及注释
如果做其他物种可以使用cellranger 的mkgtf和mkref来自己构建。

# 测试文件只要包含R1/R2即可
# 下载参考基因组及注释文件
wget https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-GRCh38-2020-A.tar.gz
# 解压
tar -zxvf refdata-gex-GRCh38-2020-A.tar.gz
# 将命令写入.sh 文件，文件内容如下
[L@loginb2 Script]$ cat CellRanger_Count.sh 
cd /Ly/exercise && \
rm -rf DEMO1 && \
cellranger count --id=DEMO1 --sample=Demo1 --fastqs=./Data/SequenceData/Demo1 --localcores=2 --localmem=20 --transcriptome=/Ly/LY_software/refdata-gex-GRCh38-2020-A --expect-cells=5000 && \
cd /Ly/exercise && \
rm -rf DEMO2 && \
cellranger count --id=DEMO2 --sample=Demo2 --fastqs=./Data/SequenceData/Demo2 --localcores=2 --localmem=20 --transcriptome=/Ly/LY_software/refdata-gex-GRCh38-2020-A --expect-cells=5000
[L@loginb2 Script]$ bash CellRanger_Count.sh

得到的结果中值得关注的文件是：
molecule_info.h5 #可用于cellranger arrg分析
filtered_feature_bc_matrix 文件夹下内容可用于下游分析
详见

[Ly@loginb2 DEMO1]$ ls outs/
analysis                       possorted_genome_bam.bam
cloupe.cloupe                  possorted_genome_bam.bam.bai
filtered_feature_bc_matrix     raw_feature_bc_matrix
filtered_feature_bc_matrix.h5  raw_feature_bc_matrix.h5
metrics_summary.csv            web_summary.html
molecule_info.h5
# filtered_feature_bc_matrix 文件夹下内容
[Ly@logina1 filtered_feature_bc_matrix]$ ls -sh
total 1.3M
 74K barcodes.tsv.gz  326K features.tsv.gz  928K matrix.mtx.gz

这三个文件分别包括什么内容？

# barcodes.tsv.gz
[Ly@logina1 filtered_feature_bc_matrix]$ zcat barcodes.tsv.gz |less -S
AAACCCAAGGAACGTC-1
AAACCCAAGGCACTAG-1
AAACCCACAAGAGGCT-1
AAACCCACACAAATCC-1
AAACCCAGTCAGGTAG-1
AAACCCAGTCGAAACG-1
AAACCCAGTTGTAAAG-1
AAACCCATCACCTTGC-1
AAACCCATCGCCTCTA-1
AAACGAAAGCGAGAAA-1
···
# features.tsv.gz
[Ly@logina1 filtered_feature_bc_matrix]$ zcat features.tsv.gz |less -S
ENSG00000243485 MIR1302-2HG     Gene Expression
ENSG00000237613 FAM138A Gene Expression
ENSG00000186092 OR4F5   Gene Expression
ENSG00000238009 AL627309.1      Gene Expression
ENSG00000239945 AL627309.3      Gene Expression
ENSG00000239906 AL627309.2      Gene Expression
ENSG00000241860 AL627309.5      Gene Expression
ENSG00000241599 AL627309.4      Gene Expression
ENSG00000286448 AP006222.2      Gene Expression
ENSG00000236601 AL732372.1      Gene Expression
···
# matrix.mtx.gz
[Ly@logina1 filtered_feature_bc_matrix]$ zcat matrix.mtx.gz |less -S
%%MatrixMarket matrix coordinate integer general
%metadata_json: {"software_version": "cellranger-6.1.1", "format_version": 2}
36601 14472 232328
5567 1 1
13754 1 1
20004 1 1
21425 1 1
26219 1 1
26563 1 2
1476 2 1
6629 2 1
10985 2 1

在下游分析中会得到一个Seurat矩阵：
barcodes.tsv.gz 是列名。
features.tsv.gz是行名。
matrix.mtx 文件是3列，第一列是行号，第二列是列号，第三列是基因表达量。
参考：表达矩阵逆转为10X的标准输出3个文件 (qq.com)
探究一下，是否真的像上述所说
下载cell Ranger运行得到的三个压缩文件这个路径下（不需要解压缩）

../Data/filtered_feature_bc_matrix/

# 创建Seurat对象
data1 <-CreateSeuratObject(Read10X("../Data/filtered_feature_bc_matrix/"))
36601x14472

矩阵大小为36601x14472和matrix.mtx.gz文件第三行的第一列、第二列一致

# matrix.mtx.gz第一列为row line ，第二列为col line
> ct=GetAssayData(object = data1, assay = "RNA", slot = "counts")
> ct[5567,1]
[1] 1
> ct[1476,2]
[1] 1
## matrix.mtx.gz第三列为表达量

matrix.mtx.gz文件压缩了所有细胞中基因表达量不为0的矩阵信息

不同样本的整合分析

一般在下游分析中整合，不用cellanger arrg

#arrg运行需要的文件
[Ly@loginb2 exercise]$ cat ./Data/Libraries.csv
sample_id,molecule_h5
Demo1,./DEMO1/outs/molecule_info.h5
Demo2,./DEMO2/outs/molecule_info.h5
# 命令
cellranger aggr --id=Demo --csv=./Data/Libraries.csv --normalize=mapped
最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 203,271评论 5赞 476
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 85,275评论 2赞 380
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 150,151评论 0赞 336
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,550评论 1赞 273
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,553评论 5赞 365
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,559评论 1赞 281
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 37,924评论 3赞 395
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,580评论 0赞 257
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 40,826评论 1赞 297
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,578评论 2赞 320
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,661评论 1赞 329
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,363评论 4赞 318
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 38,940评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,926评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,156评论 1赞 259
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 42,872评论 2赞 349
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,391评论 2赞 342