Spike in 矫正

1引言

此篇内容来自粉丝 行人 的投稿,感谢这位粉丝的分享!此外 俊粉投稿分享群 已成立,欢迎愿意分享的有志之士加入。

2投稿内容

随着高通量技术的发展，生信分析最主要的一部分内容便是组学分析。通过组学结果的分析，科研人员可以比较不同条件下的不同样本之间的差异，进而展现科学问题。然而，这种差异分析基于一个前提 “在不同实验条件下每个细胞的基因组或转录组总量是相同的”，但是这个前提很明显是错误的，即使测序前统一了上样量且测序时设定了相同的 PCR 扩增数，由于不同细胞之间的差异仍然会在文库总量上导致不可避免的“批次效应”。而目前的一些算法如单细胞 Seurat 软件取的 “Log Normalize” 的归一化和 “z-score” 标准化、bulk 转录组 DEseq2 的 Variance Stabilizing Transformation（VST） 标准化等，通过将文组标准化到一个水平来进行矫正。显然这都是无法解决这个问题的。

那么既然问题出在实验上，最好的办法还是应该回到实验上解决。早在 2014 年 ERCC（External RNA Controls Consorti） 组织为了解决这个问题提出了 spike in 的方法。Spike in 概括就是在在提取的 DNA 或 RNA 样品中加入等量的非亲缘的标准化品基因组（例如在人的建库样本中加入果蝇的基因组）简单原理如下：

图中无论是对于在全基因组范围有改变的不同样本（a），还是部分基因组区域（b）发生改变的不同样本，Spike in 矫正才能真正体现科学问题（尤其是全基因组范围的改变），而不是仅仅通过文库大小的归一化。随着二代测序技术的发展，spike in 逐渐在各种建库数据中得到了运用，并且 spike in 矫正目前被认为应该在所有的组学库中运用。（尤其对可能引起整个基因组大部分区域发生改变的实验，例如对转录因子的敲降或激活等）。

Spike in 这么重要，但生信上目前尚没有教程。这里笔者简单阐述一下关于加入 spike in 的数据分析。这里以 人的单端 RNAseq 为例，加入的是果蝇的基因组 Spike in（Chipseq 等分析在关键步奏同理，我比较喜欢用 R，所以代码全都整合到 R 了）：

# 上游
cite =”#项目目录#”
# data文件夹为存放原始下机数据的位置
files = list.files(paste0(cite,”/data”))
### 过滤以及比对
for (I in files) {
# 用来修改样本名
if(i==”AAA.gz”){name = “A-A-A”}
  system(paste0(“trim_galore -j 8 -q 30 –phred33 –length 35 –basename=”,name,” -e 0.1 –stringency 4 -o clean_data/ data/”,i))
  file.rename(from = paste0(“clean_data/”,name,”_trimmed.fq.gz”),to = paste0(“clean_data/”,name,”.fq.gz”))
  system(paste0(“fastqc -o fastqc/ -t 30 clean_data/”,name,”.fq.gz”))
  # 和人比对
  system(paste0(“bowtie2 -p 30 -x ~/gene/genomic_index/hg38 -q clean_data/”,name,”.fq.gz|samtools sort -O bam -@ 30 -l 9 > bam/”,name,”_sort.bam”))
  system(paste0(“samtools view -h -F 4 -q 30 -@ 30 bam/”,name,”_sort.bam|samtools sort -O bam -@ 30 -l 9 > bam/”,name,”.bam”))
  # 和果蝇spike in比对
  system(paste0(“bowtie2 -p 30 -x /reference/dm6/bowtie2/genome -q clean_data/”,name,”.fq.gz |samtools sort -O bam -@ 30 -l 9 > bam_sbink/”,name,”.bam”))
  system(paste0(“samtools index -@ 30 -b bam/”,name,”.bam”))
  }

system(paste0(“multiqc fastqc/ -f -n clean -o ./multiqc”))
### 前面的分析除需要增加和果蝇的比对外，和普通的RNAseq无异

### 生成bw文件
if(!dir.exists(“bw”)){dir.create(“bw”)}
files = list.files(paste0(cite,”/bam”),pattern = “.bam”) %>% grep(“sort”,.,invert = T,value = T)  %>% grep(“.bai”,.,invert = T,value = T)

setwd(paste0(cite,”/bam”))
# 通过spike in计算标准化因子
system(paste0(“multiBamSummary bins –bam “,paste(files,collapse = “ “),” –scalingFactors scalingfactors.txt -p 30”))

scale = read.table(“scalingfactors.txt”,row.names = 1,header = 1)

setwd(cite)
for (I in rownames(scale)) {
#这些文件就可以扔到IGV去查看了
  system(paste0(“bamCoverage -p 30 –normalizeUsing RPKM –scaleFactor “,scale[I,”scalingFactor”],” -b bam/”,I,” -o bw/”,sub(“.bam”,””,i),”.bw”))
  }

# 下游
## 此外，除了上游出bw文件需要用到sbike in外，下游的DESeq2包中计算标准化因子需要用到，这里主要展示构建dds的代码其他分析和普通RNAseq一致：
dds = DESeqDataSetFromMatrix(countData = as.matrix(countdata), colData=database, design= ~ condition)
spike = read.table(paste0(cite,"/bam/scalingfactors.txt"),sep="t",header=T,row.names = 1)
rownames(spike) = sub("\..*$","",rownames(spike))
# 先让他自己算
dds = estimateSizeFactors(dds)
DESeq2::sizeFactors(dds) = as.numeric(1/spike[names(DESeq2::
# 咱不用它算的，用咱们自己的
sizeFactors(dds)),"scalingFactor"])
dds <- DESeq(dds)
mcols(dds)$basepairs = as.numeric(row_data$Length)
dds <- dds[rowMaxs(fpkm(dds))>1,]

总的来说 spike in 对我们的科研还是很有助益的，现在一些顶尖的文章中例如 Chipseq、TTseq 等都会使用该方法来提高实验的准确性。而生信分析也不难，主要是抓住几个标准化的点。

欢迎来投稿!

欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 (微信交流群需收取 20 元入群费用,一旦交费,拒不退还!(防止骗子和便于管理)) 。QQ 群可免费加入, 记得进群按格式修改备注哦。