「BioNano系列」那些Bionano未覆寫的區域是什麼？

在 「Bionano系列」光學圖譜混合組裝應該怎麼做？

這篇文章中，我展示了下面這張圖。

和之前的圖不同的是，我加了幾個箭頭，這些箭頭所指向的區域的特征就是，這些區域并未被Bionano所覆寫。如果不去思考這些區域到底是什麼，直接進行混合組裝，那麼這其實對最後結果的不負責任。因為這完全可能是組裝軟體沒有正确的處理錯誤的overlap，将不應該連接配接的序列連接配接在一起（盡管這個機率不高）。

我的直覺猜測就是，這些區域應該是重複序列區域。畢竟Bionano标記技術依賴于酶識别特定位點進行酶切加上熒光标記，重複序列要麼會因酶切密度太高，相機的分辨率達不到，而識别失敗，要麼是酶切位點過少，信号太弱。

那麼我應該如何驗證這個猜想？通過幾天的文獻翻閱和嘗試，我用重複序列數量和基因數量的相對比值進行衡量。

指令行的代碼如下(沒有考慮檔案的相對位置)

# 利用拟南芥的原本CDS進行注釋
gmap_build -D index -d R05C0144 ../R05C0144.fa &
gmap -t 20 -D index -d R05C0144 -f gff3_gene ../Athaliana_cds.fa > cds_gene.gff3 2> log.txt &
# 重複序列注釋
RepeatMasker -e ncbi -species arabidopsis -pa 30 -gff -dir . ../R05C0144.fa &
# GFF轉成BED
awk 'BEGIN{OFS="\t"} {print $1,$4,$5}' ../repeat_annotation/R05C0144.fa.out.gff > repeat.bed
grep -w 'gene' ../gene_annotation/cds_gene.gff3| awk 'BEGIN{OFS="\t"} {print $1,$4,$5}' | bedtools sort -i - > gene.bed
# 統計
bedtools makewindows -w 100000 -g ../R05C0144.txt > windows_100k.bed
bedtools coverage -a windows_100k.bed -b repeat.bed > repeat_stat.bed
bedtools coverage -a windows_100k.bed -b gene.bed > gene_stat.bed           

R代碼如下

gene_df <- read.table("R05C0144/feature_stat/gene_stat.bed",
                        sep="\t", stringsAsFactors = F)
repeat_df <- read.table("R05C0144/feature_stat/repeat_stat.bed",
                      sep="\t", stringsAsFactors = F)

options(scipen=999) 
contig <- "contig2"

repeat_ctg <- repeat_df[repeat_df$V1 == contig,]
gene_ctg <- gene_df[gene_df$V1 == contig,]

combine_df <- data.frame(pos=(repeat_ctg$V2 + repeat_ctg$V3) / 2,
                         repeat_num=repeat_ctg$V4,
                         gene_num=gene_ctg$V4)
combine_df$total = combine_df$repeat_num + combine_df$gene_num

combine_df$gene_ratio <- combine_df$gene_num / combine_df$total * 100

combine_df$repeat_ratio <- combine_df$repeat_num / combine_df$total * 100


plot(combine_df$pos, combine_df$gene_ratio, 
     type="l", 
     ylim=c(0,100),
     xlab="position",
     ylab="percent",
     col="blue")
lines(combine_df$pos, combine_df$repeat_ratio, col="red")
abline(v=7.85*1e6)           

我檢查了一些區間，的确是重複序列比例高于基因比例，當然還有一些區間不是。說明重複序列并不是光學圖譜未覆寫的主要原因。

當然對于拟南芥這種有着高品質基因組的物種而言，我們還可以進行共線性分析。不過對于這些N50在4M左右，而且低雜合的基因組，其實都不需要太操心這種錯誤。

我這裡也就驗證了一種可能性，後續還得檢查了一下其他原因，說不定僅僅是光學圖譜的深度不夠而已。