Sam 文件 flag 研究

原创 JunJunLab 老俊俊的生信笔记 2022-08-17

收录于合集 #Ribo-seq 39个

要怎么在这薄情的世界里,深情地活着?

1引言

为什么会与人研究这玩意呢? 因为在我后续的分析中,忽略了这个 flag 所带来的影响,导致我分析的结果完全不一样。

之前比对到 转录组序列 上的 sam 文件用来做质控,检查 Riboseq 的数据质量,结果如下:

我后面用同样的 fastq 文件 比对到基因组, 写了类似的代码来质控,结果如下:

是不是完全不一样了,于是就开启了后面的探索过程。

2sam 文件探索

我们先看看比对到基因组序列上的 sam 文件的 flag 组成:

$ less -S WT_Ribo2.sam|grep -v '@'|awk '{print $2}'|sort|uniq -c
1034988 0
1595628 16
2124062 4

根据 sam 文件 flag 值定义:

flag0: 单端测序,且该 reads 与序列相同。
flag16: 与参考序列反向互补。
flag4: 未比对上的 reads。

由以上结果可以看到,对于比对上的 reads 来说,有快一半的序列是 flag0,有一半多的序列为 flag16。于是我分别提取出 flag0 和 flag16 的序列分别进行质控:

# 提取flag 0的序列
$ mtools view -F 20 -O SAM WT_Ribo2.sam > WT_Ribo2_0.sam

# 提取flag 16 的序列
$ samtools view -f 16 -O SAM WT_Ribo2.sam > WT_Ribo2_16.sam

分析跳过,拿分析结果绘图:

library(tidyverse)
library(ggplot2)
library(cowplot)

# load data
file <- list.files(pattern = 'error.txt')

map_df(1:length(file),function(x){
  tmp = read.table(file[x])
  colnames(tmp) <- c('length','framest','relst','framesp','relsp','counts')
  tmp$sample <- sapply(strsplit(file[x],split = '\\.'),'[',2)
  return(tmp)
}) -> dfqc

#################################################
# frame
#################################################

frame <- dfqc %>% group_by(sample,length,framest) %>% summarise(num = sum(counts))

# plot
ggplot(frame,
  aes(x = factor(length),y = num/1000)) +
  geom_col(aes(fill = factor(framest)),position = position_dodge2()) +
  scale_fill_brewer(palette = 'Greens',name = '',direction = -1) +
  theme_bw(base_size = 16) +
  theme(aspect.ratio = 0.8,
        strip.background = element_rect(colour = NA),
        legend.position = 'top',
        legend.background = element_blank(),
        panel.grid = element_blank()) +
  xlab('Read length') + ylab('Reads numbers(k)') +
  facet_wrap(~sample,scales = 'free',ncol = 2)

可以看到只有 flag16 的质控结果和比对到转录组上的结果才一致。

为了进一步探究原因我们查看比对到转录组上的 sam 文件的 flag 组成:

$ less -S WT_Ribo2.sam|grep -v '@'|awk '{print $2}'|sort|uniq -c
  13697 0
1677169 16
3063812 4

可以看到 flag16 的占百分之 90 多,所以整体效果会好一些,因为 flag0 的扰动小了。

flag0 到底是个啥?

我们知道 基因组序列和处在正链的基因的序列是一致 , 负链的则是反向互补。

Ribo-seq 的过程是捕获核糖体保护的 mRNA 片段,然后对这些片段进行 反转录进行测序, 测序的片段如果比对到基因组上面,那么 flag16(即与基因组反向互补)的序列应该是 正链的基因, 而 flag0(与基因组序列相同)的序列应该是 负链的基因。

我们对 flag0 和 flag16 的 sam 转为 bw 文件在 IGV 里进行查看:

$ samtools sort -@ 10 WT_Ribo2.sam -o WT_Ribo2.bam
$ samtools index WT_Ribo2.bam

# flag16 to bigwig
$ samtools view -f 16 WT_Ribo2.bam -O BAM > flag16.bam
$ samtools index flag16.bam
$ bamCoverage --normalizeUsing BPM -p 10 -b flag16.bam -o flag16.bw

# flag0 to bigwig
$ samtools view -F 20 WT_Ribo2.bam  -O BAM > flag0.bam
$ samtools index flag0.bam
$ bamCoverage --normalizeUsing BPM -p 10 -b flag0.bam -o flag0.bw

由以上结果可以清楚的看到和我说的结论是一致的。

但为什么 flag0 的质控结果不一样呢?

这里只能锁定是代码的原因了,我写的质控代码并没有把 flag0 和 flag16 考虑进去,计算 frame 都是拿 reads 的 5'end 去计算离 start codon/stop codon 的距离来分析的。

对于正链基因(flag16)的代码应该是正确的,我们只要分析负链基因(flag0)的算法该如何改正就行了。根据 Ribo-seq 原理及比对的 flag0,对于负链基因,sam 文件的比对的位置虽然是 5'end 的位置,但是却是远离负链基因的翻译反向那一端的,示意图如下:

所以计算 frame 我们应该使用 3'end 来计算:

3修改代码

理解了上面的东西,我们只需要将 flag 考虑进去,对于 flag0 的则用 3'end 来计算,而 flag16 的则用 5'end 来分析。

代码修改部分:(julia语言)

# tags
refname,align_pos,read_length = SAM.refname(record),SAM.position(record),SAM.seqlength(record)

# read flag tag
flag = SAM.flag(record)

# flag16(+) use 5'end as alignpos and flag0(-) use 3'end as alignpos
if flag == 16
    # flag 16 reads from + stand
    # read key
    readKey = join([refname,align_pos],"|")
elseif flag == 0
    # flag 0 reads from - stand
    end3Pos = align_pos + read_length - 1
    # read key
    readKey = join([refname,end3Pos],"|")
end