实战演练 理论知识学再好,能付诸实践灵活运用才行,所以我们常强调知行合一,实践出真知。实战演练这个栏目就是带大家从头到尾完整复现单细胞文献分析流程。好了,干货多,屁话少,我们来看实战流程。 希望大家能有所收获!
目的:将每个flowcell 的Illumina sequencer's base call files (BCLs)转为fastq文件 特色: 它借鉴了Illumina出品的
它的示意流程: 两种使用方式: $ cellranger mkfastq --id=bcl \ --run=/path/to/bcl \ --samplesheet=samplesheet-1.2.0.csv # 第二种 $ cellranger mkfastq --id=bcl \ --run=/path/to/bcl \ --csv=simple-1.2.0.csv # 其中id指定输出目录的名称,run指的是下机的原始BCL文件目录 # 重要的就是测序lane、样本名称、index等信息 samplesheet.csv文件就是illumina常规使用的,类似下面这种。它除了需要指定各种ID、name之外,还要根据不同的试剂盒版本调整
还有一种10X定制的简单化的csv文件,例如: 1,test_sample,SI-GA-A3 # 其中第一列指定lane ID,第二列是样本名称,第三列是index名称 使用简化版的这个文件,可以识别使用的试剂盒版本,然后自行调整reads的长度信息
目录结构如下: - tiny-bcl/outs/fastq_path/bcl/- Sample1 - Sample1_S1_L001_I1_001.fastq.gz - Sample1_S1_L001_R1_001.fastq.gz - Sample1_S1_L001_R2_001.fastq.gz 自己分析的数据也要改成这种结构存放,方便后续分析
第一种情况:利用 # 1.所有mkfastq生成的样本 --fastqs=MKFASTQ_ID/outs/fastq_path # 2. 多个flowcell生成的所有样本 --fastqs=MKFASTQ_ID/outs/fastq_path1,MKFASTQ_ID/outs/fastq_path2 # 3.所有bcl2fastq 生成的样本 --fastqs=/PATH/TO/bcl2fastq_output # 4. 所有lanes上的test_sample1样本 --fastqs=MKFASTQ_ID/outs/fastq_path \ --sample=test_sample1 # 5. lane1上的test_sample1样本 --fastqs=MKFASTQ_ID/outs/fastq_path \ --sample=test_sample1 \ --lanes=1 # 6. 将test_sample1和test_sample2各自进行操作 --fastqs=MKFASTQ_ID/outs/fastq_path \ --sample=test_sample1,test_sample2 其实从上面的各种设置也能看出来,一开始的样本命名规则是非常重要的 第二种情况:也是利用 --fastqs=MKFASTQ_ID/outs/fastq_path \ --sample=SI-GA-A1_1,SI-GA-A1_2,SI-GA-A1_3,SI-GA-A1_4 # 2. 只处理SI-GA-A1样本的第一个index样本 --fastqs=MKFASTQ_ID/outs/fastq_path \ --sample=SI-GA-A1_1 第三种情况:也是利用 --fastqs=MKFASTQ_ID/outs/fastq_path # 2. bcl2fastq得到的所有样本 --fastqs=/PATH/TO/bcl2fastq_output # 3. test_sample样本的所有lanes --fastqs=MKFASTQ_ID/outs/fastq_path \ --sample=test_sample # 4. test_sample样本的lane1 --fastqs=MKFASTQ_ID/outs/fastq_path \ --sample=test_sample \ --lanes=1 第四种情况:使用 --fastqs=/PATH/TO/PROJECT_FOLDER # 2. 处理Mysample样本的所有lanes的数据 --fastqs=/PATH/TO/PROJECT_FOLDER \ --sample=MySample # 3. 只处理Mysample样本的lane1数据 --fastqs=/PATH/TO/PROJECT_FOLDER \ --sample=MySample \ --lanes=1 第五种情况:fastq命名方式变了,类似于这样: 它一般是从 --fastqs=/PATH/TO/PROJECT_FOLDER # 2.所有SI-GA-A1样本 --fastqs=/PATH/TO/PROJECT_FOLDER \ --indices=SI-GA-A1 # 3.所有SI-GA-A1样本的lane1数据 --fastqs=/PATH/TO/PROJECT_FOLDER \ --indices=SI-GA-A1 \ --lanes=1 # 4.利用oligo --fastqs=/PATH/TO/PROJECT_FOLDER \ --indices=AACCGTAA,CTAAACGG,GGTTTACT,TCGGCGTC 第六种情况:数据命名与上面完全不同,因此需要自己重命名,方式就是 # 这个在单细胞实战(二)中介绍过[Sample Name]_S1_L00[Lane Number]_[Read Type]_001.fastq.gz # 其中Read Type # I1: Sample index read (optional) # R1: Read 1 # R2: Read 2 分析时就可以直接调用了 --fastqs=/PATH/TO/PROJECT_FOLDER # 2. 某个样本的所有lanes数据 --fastqs=/PATH/TO/PROJECT_FOLDER \ --sample=SAMPLENAME # 3. 某个样本的某个lane --sample=SAMPLENAME \ --fastqs=/PATH/TO/PROJECT_FOLDER \ --lanes=1 这个过程是最重要的,它完成细胞与基因的定量,它将比对、质控、定量都包装了起来,内部流程很多,但使用很简单 先学会使用每个版本要求的参数是不同的,尤其是V2与V3版本存在较大差异,这里先对V2进行了解 基本上自己需要输入的参数是: # 这是示例,不是真实数据 #cellranger count --id=sample345 \ --transcriptome=/opt/refdata-cellranger-GRCh38-1.2.0 \ --fastqs=/home/scRNA/runs/HAWT7ADXX/outs/fastq_path \ --sample=mysample \ --expect-cells=1000 \ --nosecondary # id指定输出文件存放目录名 # transcriptome指定与CellRanger兼容的参考基因组 # fastqs指定mkfastq或者自定义的测序文件 # sample要和fastq文件的前缀中的sample保持一致,作为软件识别的标志 # expect-cells指定复现的细胞数量,这个要和实验设计结合起来 # nosecondary 只获得表达矩阵,不进行后续的降维、聚类和可视化分析(因为后期会自行用R包去做) 它的输出文件有很多Outputs:- Run summary HTML: /opt/sample345/outs/web_summary.html - Run summary CSV: /opt/sample345/outs/metrics_summary.csv - BAM: /opt/sample345/outs/possorted_genome_bam.bam - BAM index: /opt/sample345/outs/possorted_genome_bam.bam.bai - Filtered gene-barcode matrices MEX: /opt/sample345/outs/filtered_gene_bc_matrices - Filtered gene-barcode matrices HDF5: /opt/sample345/outs/filtered_gene_bc_matrices_h5.h5 - Unfiltered gene-barcode matrices MEX: /opt/sample345/outs/raw_gene_bc_matrices - Unfiltered gene-barcode matrices HDF5: /opt/sample345/outs/raw_gene_bc_matrices_h5.h5 - Secondary analysis output CSV: /opt/sample345/outs/analysis - Per-molecule read information: /opt/sample345/outs/molecule_info.h5 - Loupe Cell Browser file: /opt/sample345/outs/cloupe.cloupe Pipestance completed successfully! 从上到下依次来看:
一些内置软件和算法基因组比对—是否在外显子?利用了 STAR比对工具,这款比对工具比对速度快,灵敏度高,是ENCODE、GATK推荐使用的工具,允许基因的可变剪切。比对完之后,利用GTF文件将reads溯源回外显子区、内含子区、基因间区:如果一条read的50%以上与外显子有交集,那么就认为它在外显区;如果不在外显子区,与内含子有交集,那么就认为它在内含子区;与外显子、内含子都没有交集,那么就认为在基因间区 MAPQ 辅助判断—在外显子的正确率有多少?如果reads比对到了一个外显子区,同时也比对到了1个或多个的非外显子区,更相信它在外显子区,然后看MAPQ值,值越大越可信,如果MAPQ的值为255的话,那么就可以非常确定它比对到了外显子区
转录组比对—是否特异比对?如果上面得到的外显子区域reads同时比对上有注释转录本上的外显子,并且在同一条链上,那么认为这个reads也比对到了转录组;如果只比对到单个基因的注释信息,那么认为它是特异比对到转录组的(uniquely /confidently mapped ),这样的reads才会拿来做接下来的UMI 计数 重点和难点在于自主构建参考信息Cell Ranger为比对和定量提供了参考基因组及注释 pre-built human (hg19, GRCh38), mouse (mm10), and ercc92 reference packages 但是很多时候,我们需要根据自己的需要,自定义一套参考信息,但需要注意以下问题:
首先利用mkgtf过滤GTF文件先从 ENSEMBL或UCSC上下载,然后使用 # 其中键值对可以指定多个,比如 $ cellranger mkgtf Homo_sapiens.GRCh38.ensembl.gtf Homo_sapiens.GRCh38.ensembl.filtered.gtf \ --attribute=gene_biotype:protein_coding \ --attribute=gene_biotype:lincRNA \ --attribute=gene_biotype:antisense \ --attribute=gene_biotype:IG_LV_gene \ --attribute=gene_biotype:IG_V_gene \ --attribute=gene_biotype:IG_V_pseudogene \ --attribute=gene_biotype:IG_D_gene \ --attribute=gene_biotype:IG_J_gene \ --attribute=gene_biotype:IG_J_pseudogene \ --attribute=gene_biotype:IG_C_gene \ --attribute=gene_biotype:IG_C_pseudogene \ --attribute=gene_biotype:TR_V_gene \ --attribute=gene_biotype:TR_V_pseudogene \ --attribute=gene_biotype:TR_D_gene \ --attribute=gene_biotype:TR_J_gene \ --attribute=gene_biotype:TR_J_pseudogene \ --attribute=gene_biotype:TR_C_gene # 这样得到的Homo_sapiens.GRCh38.ensembl.filtered.gtf结果中就不包含gene_biotype:pseudogene这部分 然后利用mkref构建参考索引# 基本使用(单个物种)cellranger mkref --genome=hg19 --fasta=hg19.fa --genes=hg19-filtered-ensembl.gtf # 可以使用--nthreads使用多线程加速 # 得到的输出结果(保存在--genome名称的目录中) ls hg19 fasta/ genes/ pickle/ reference.json star/ # 如果对于多个物种组合(本文的数据其实就应该这样组合起来) cellranger mkref --genome=hg19 --fasta=hg19.fa --genes=hg19-filtered-ensembl.gtf \ --genome=mm10 --fasta=mm10.fa --genes=mm10-filtered-ensembl.gtf # 得到的结果 ls hg19_and_mm10 fasta/ genes/ pickle/ reference.json star/ 如果要增加基因信息参考链接:https://kb./hc/en-us/articles/115003327112-How-can-we-add-genes-to-a-reference-package-for-Cell-Ranger- 第一步,在 第二步,在 # 第一列:Chromosome 指定基因组上染色体或contig位置 # 第二列:Source 这个用处不大 # 第三列:Feature CellRanger软件只取exon的部分 # 第四列:Start 起始位点(1-based) # 第五列:End 终止位点(1-based) # 第六列:Score 这个用处不大,建议用"."表示 # 第七列:Strand feature信息在基因组的+或-链 # 第八列:Frame 用处不大,建议“.” # 第九列:分号分隔的键值对,重点是transcript_id 和gene_id。gene_name可选 例如: mylocus annotation exon 100 200 . + . gene_id "mygene"; transcript_id "mygene"; 第三步,使用 P.S. 最后得到的参考信息(包括参考基因组、注释信息)文件结构如下:# 这是官网下载的hg38数据refdata-cellranger-GRCh38-1.2.0 |-- [ 222] README.BEFORE.MODIFYING |-- [4.0K] fasta | `-- [2.9G] genome.fa |-- [4.0K] genes | `-- [1.3G] genes.gtf |-- [4.0K] pickle | `-- [ 58M] genes.pickle |-- [ 424] reference.json |-- [4.0K] star | |-- [3.0G] Genome | |-- [8.0G] SA | |-- [1.5G] SAindex | |-- [1.2K] chrLength.txt | |-- [1.9K] chrName.txt | |-- [3.0K] chrNameLength.txt | |-- [2.1K] chrStart.txt | |-- [ 37M] exonGeTrInfo.tab | |-- [ 15M] exonInfo.tab | |-- [526K] geneInfo.tab | |-- [ 909] genomeParameters.txt | |-- [9.1M] sjdbInfo.txt | |-- [7.1M] sjdbList.fromGTF.out.tab | |-- [7.1M] sjdbList.out.tab | `-- [9.4M] transcriptInfo.tab `-- [ 6] version 4 directories, 21 files 当处理多个生物学样本或者一个样本存在多个重复/文库时,最好的操作就是先分别对每个文库进行单独的count定量,然后将定量结果利用 第一步 得到count结果例如现在分别进行3个定量流程 $ cellranger count --id=LV123 ...... wait for pipeline to finish ... $ cellranger count --id=LB456 ... ... wait for pipeline to finish ... $ cellranger count --id=LP789 ... ... wait for pipeline to finish ... 第二步 构建Aggregation CSV就像这样: # AGG123_libraries.csvlibrary_id,molecule_h5 LV123,/opt/runs/LV123/outs/molecule_info.h5 LB456,/opt/runs/LB456/outs/molecule_info.h5 LP789,/opt/runs/LP789/outs/molecule_info.h5 # 其中 # molecule_h5:文件molecule_info.h5 file的路径 第三步 运行aggrcellranger aggr --id=AGG123 \--csv=AGG123_libraries.csv \ --normalize=mapped # 结果输出到AGG123这个目录中 至于最后的 reanalyze ,这个属于定制化分析,这里暂时不做探讨,日后待标准化流程构建起来,再补充这一部分 |
|