scRNA-seq数据处理—文件格式小结

书籍翻译

好的书籍是人类进步的阶梯，但有些人却找不到优秀的阶梯，为此我们开设了书籍翻译这个栏目，作为你学习之路的指路明灯；分享国内外优秀书籍，弘扬分享精神，做一个知识的传播者。

希望大家能有所收获！

目录

正

文

处理原始scRNA-seq数据

3.3

文件格式

3.3.1 FastQC

      FastQ是您将遇到的最原始形式的scRNASeq数据。所有scRNASeq方案都使用配对末端测序进行测序。Barcode序列可以在一个或两个reads中发生，这取决于所采用的protocol 。然而，使用独特分子标识符（UMI）的protocol 通常包含一个带有细胞和UMI barcode 和 adapters 但没有任何转录序列的read。因此，尽管实际上是成对末端测序，但reads将被比对为好像它们是单端测序的。

FastQ文件的格式如下：

>ReadID
READ SEQUENCE
+
SEQUENCING QUALITY SCORES

3.3.2 BAM

      BAM文件格式以标准且有效的方式存储比对过的reads。人类可读的版本称为SAM文件，而BAM文件是高度压缩的版本。BAM / SAM文件包含标题。标题通常包括有关样品制备，测序和比对的信息; 和每个read的每个比对的制表符分隔行。

alignment行使用具有以下列的标准格式：

1. QNAME：read名称（通常包括UMI条形码）

2. FLAG：数字标记表示比对的“类型”，链接：所有可能的“类型”的解释

3. RNAME：参考序列名称（即染色体读数被比对到了什么序列上）

4. POS：最左边的比对位置

5. MAPQ：比对质量

6. CIGAR：read的匹配/不匹配部分的字符串（可能包括soft-clipping）

7. RNEXT：配对/下个read的参考名称

8. PNEXT：配对/下个read的POS

9. TLEN：模板长度（read被比对到的参考区域的长度）

10. SEQ：read序列

11. QUAL：read质量

可以使用samtools将BAM / SAM文件转换为其他格式：

samtools view -S -b file.sam > file.bam
samtools view -h file.bam > file.sam

一些测序设备会自动将您的read比对到标准基因组，并提供BAM或CRAM格式的文件。通常它们不会在基因组中包含ERCC序列，因此在BAM / CRAM文件中不会比对ERCC read。要量化ERCC（或任何其他遗传变化），或者如果您只想使用不同于通用pipeline中的任何比对算法（通常是过时的算法），那么您需要将BAM / CRAM文件转换回FastQs：

可以使用bedtools将BAM文件转换为FastQ。为了确保多比对reads的单个拷贝首先按read名称排序，并使用samtools删除次级比对。Picard也包含了一种将BAM转换为FastQ文件的方法。

# sort reads by name
samtools sort -n original.bam -o sorted_by_name.bam
# remove secondary alignments
samtools view -b -F 256 sorted_by_name.bam -o primary_alignment_only.bam
# convert to fastq
bedtools bamtofastq -i primary_alignment_only.bam -fq read1.fq -fq2 read2.fq

3.3.3 CRAM

      CRAM文件类似于BAM文件，只有它们包含用于header中比对到的参考基因组的信息。这允许在每个read里的与参考基因组相同的碱基被进一步压缩。与BAM相比，CRAM还支持一些有损（lossy）数据压缩的方法以进一步优化存储。CRAM主要由Sanger / EBI测序设备使用。

CRAM和BAM文件可以使用最新版本的samtools（> = v1.0）进行转换。但是，这种转换可能需要将参考基因组下载到缓存（cache）中。或者，您可以从CRAM文件的header中的元数据（metadata）预先下载正确的参考基因组，或者通过与生成CRAM的人交谈，并使用'-T'指定该文件，因此我们建议在执行此操作之前设置特定的缓存位置：

export REF_CACHE=/path_to/cache_directory_for_reference_genome
samtools view -b -h -T reference_genome.fasta file.cram -o file.bam
samtools view -C -h -T reference_genome.fasta file.bam -o file.cram

3.3.4 手动检查文件

      有时，手动检查文件可能很有用，例如检查文件来自正确样本的header中的元数据。'less'和'more'可用于检查命令行中的任何文本文件。通过使用“|”将samtools视图的输出到这些命令中，而不必保存每个文件的多个副本。

less file.txt
more file.txt
# counts the number of lines in file.txt
wc -l file.txt
samtools view -h file.[cram/bam] | more
# counts the number of lines in the samtools output
samtools view -h file.[cram/bam] | wc -l

练习

您已经获得了一个小的cram文件：EXAMPLE.cram

任务1：此文件是如何比对出来的？使用了什么软件？使用了什么基因组？（提示：检查header）

任务2：多少reads被比对了/没被比对？总reads数是多少？有多少次级比对？（提示：使用FLAG）

任务3：将CRAM转换为两个Fastq文件。每个read都得到一份拷贝吗？（将这些文件命名为“10cells_read1.fastq”“10cells_read2.fastq”）

如果您遇到困难，可以通过输入运行命令naked来显示每个软件的帮助信息 - 例如'samtools view'，'bedtools'

3.3.5 基因组（FASTA GTF）

      要比对您的reads，您还需要参考基因组，在许多情况下还需要基因组注释文件（采用GTF或GFF格式）。这些可以从任意的主要基因组学数据库下载：Ensembl，NCBI或UCSC Genome Browser。

GTF文件包含基因，转录本和外显子的注释。它们一定包含：（1）seqname：chromosome / scaffold（2）source：这个注释来自哪里（3）feature：这是什么类型的特征？（例如基因，转录本，外显子）（4）start：开始位置（bp）（5）end：结束位置（bp）（6）score：数字（7）strand：+（前进）或 - （反向）（ 8）frame：CDS指示哪个碱基是第一个密码子的第一个碱基（0 =第一个碱基，1 =第二个碱基等等）。（9）attribute：以分号分隔的标签值对的额外信息对的列表（例如姓名/身份证，生物类型）

空字段标有“。”。

根据我们的经验，Ensembl是最容易使用的，并且具有最大的注释集。NCBI往往更严格，仅包括高置信度基因注释。而UCSC包含多个使用不同标准的基因组注释。

如果您的实验系统包含非标准序列，则必须将这些序列添加到基因组fasta和gtf中以量化它们的表达。最常见的是，这是针对ERCC加标进行的，尽管必须对CRISPR相关序列或其他过表达/报告构建体进行相同的操作。

为了获得最大的有效性/灵活性，我们建议为所有非标准序列创建完整和详细的entries。

没有标准化的方法来做到这一点。以下是我们的自定义perl脚本，用于为ERCC创建一个gtf和fasta文件，可以将其附加到基因组中。当/如果要量化内含子reads时，您可能还需要更改gtf文件以处理内含子中的重复元素。任何脚本语言甚至“awk”或一些文本编辑器都可以用来相对有效地完成这项任务，但它们超出了本课程的范围。

# Converts the Annotation file from 
# https://www.thermofisher.com/order/catalog/product/4456740 to 
# gtf and fasta files that can be added to existing genome fasta & gtf files.

my @FASTAlines = ();
my @GTFlines = ();
open (my $ifh, "ERCC_Controls_Annotation.txt") or die $!;
<$ifh>; #header
while (<$ifh>) {
  # Do all the important stuff
  chomp;
  my @record = split(/\t/);
  my $sequence = $record[4];
  $sequence =~ s/\s+//g; # get rid of any preceeding/tailing white space
  $sequence = $sequence."NNNN";
  my $name = $record[0];
  my $genbank = $record[1];
  push(@FASTAlines, ">$name\n$sequence\n");
# is GTF 1 indexed or 0 indexed? -> it is 1 indexed
# + or - strand?
  push(@GTFlines, "$name\tERCC\tgene\t1\t".(length($sequence)-2)."\t.\t+\t.\tgene_id \"$name-$genbank\"; transcript_id \"$name-$genbank\"; exon_number \"1\"; gene_name \"ERCC $name-$genbank\"\n");
  push(@GTFlines, "$name\tERCC\ttranscript\t1\t".(length($sequence)-2)."\t.\t+\t.\tgene_id \"$name-$genbank\"; transcript_id \"$name-$genbank\"; exon_number \"1\"; gene_name \"ERCC $name-$genbank\"\n");
  push(@GTFlines, "$name\tERCC\texon\t1\t".(length($sequence)-2)."\t.\t+\t.\tgene_id \"$name-$genbank\"; transcript_id \"$name-$genbank\"; exon_number \"1\"; gene_name \"ERCC $name-$genbank\"\n");
} close($ifh);

# Write output
open(my $ofh, ">", "ERCC_Controls.fa") or die $!;
foreach my $line (@FASTAlines) {
  print $ofh $line;
} close ($ofh);

open($ofh, ">", "ERCC_Controls.gtf") or die $!;
foreach my $line (@GTFlines) {
  print $ofh $line;
} close ($ofh);

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

• 全国巡讲第6、7站-郑州和西安（生信技能树爆款入门课）（火热报名中）
  

单细胞天地欢迎你

单细胞天地

生信技能树