CheatSheet

首先先放上 0-based 和 1-based 的cheatsheet

關(guān)于USCS瀏覽器背后的數(shù)據(jù)集，可見Database/browser start coordinates differ by 1 base

關(guān)于SAF的格式，網(wǎng)上好像沒有找到確切的答案。只有 biostar 這里提了下，認為是0-based

IGV的話，至少從fa的坐標來看，是1-based的

對于BigWig格式而言，有兩種類型，具體見1-start coordinate system in use for variableStep and fixedStep

• 如果其是從bedgraph創(chuàng)建而來的，那么其就是0-based。
• 如果是從wiggle創(chuàng)建而來，其就是1-based。

cheatsheet的整合綜合了

• The UCSC Genome Browser Coordinate Counting Systems
• Tutorial: Cheat Sheet For One-Based Vs Zero-Based Coordinate Systems
• Bioinformatics Data Skill-p268的表格

1-based與0-based

事實上，生信的文件系統(tǒng)中存在著兩種坐標系統(tǒng)，1-based 和 0-based。這兩種坐標系統(tǒng)各有利弊，產(chǎn)生的原因也很復雜，其中還關(guān)系到不同編程語言的坐標格式。0-based和1-based的主要區(qū)別和優(yōu)勢有

0-based的話，其計算區(qū)間更加的方便，比如一個區(qū)間是[4,10)，那么其就是6個堿基長度。而對應的，1-based，就是[5,10]了，要 10-5+1 = 6 。

0-based事實上還可以代表 寬度為0 的區(qū)間，這一點在你定義核酸酶切割位點的時候是比較有用的，比如[23, 23)。畢竟你切得是兩個堿基之間，而不是某個特定堿基。當然 1-based 的話，有時候會用[24,23] 來代表……

Py（應該還有C）的話，對于一個字符串的選擇是從0開始的，其在區(qū)間選擇的時候，是半閉半開區(qū)間，比如下面

>>> "CTTACTTCGAAGGCTG"[1:5]
'TTAC'

而R在區(qū)間選擇的時候，對于一個字符串的選擇是從1開始的，用到的是閉區(qū)間

> substr("CTTACTTCGAAGGCTG", 2, 5)
[1] "TTAC"

這部分主要來自Bioinformatics Data Skill p265-266

事實上，對于 0-based 和 1-based的解釋，我更喜歡 Tutorial: Cheat Sheet For One-Based Vs Zero-Based Coordinate Systems 這篇文章里面的。

對于像下面從某條染色體起始的7個堿基

image

• 1-based：數(shù)字是直接代表一個堿基的
• 0-based：在兩個數(shù)字之間才代表一個堿基

image

• 1-based：單個堿基，單個多態(tài)性位點，堿基區(qū)間都是與數(shù)字對應的
• 0-based：單個堿基，單個多態(tài)性位點，堿基區(qū)間是由兩個數(shù)字之間對應的

image

• 1-based：
  • 缺失對應的是這個缺失堿基的數(shù)字坐標，這里是5號堿基T的缺失，我們之前使用[5, 5]代表的T，那么坐標就是[5, 5]
  • 插入對應的是兩個堿基數(shù)字坐標中間，這里是3號C和4號G之間插入了TTA。我們之前使用[3, 4]代表的CG，那么插入坐標就是[3,4]
• 0-based：
  • 缺失對應的是這個缺失堿基兩側(cè)的數(shù)字坐標，比如我們之前是用[4,5)代表的T，那么這里就是4-5
  • 插入對應的是這個插入發(fā)生的數(shù)字坐標，這里插入發(fā)生在C和G之間。我們用[2,3)代表C，用[3,4)代表了G，那么[3,3)就代表了CG中間那個位置。

其實 0-based 和 1-based 一般還不太要緊，但對于 VCF 這種錙銖必較的位點位點，了解 0 和 1就很重要了。

實際上，0 和 1之間的轉(zhuǎn)換還是比較方便的，根據(jù) The UCSC Genome Browser Coordinate Counting Systems 里面說

• 0-start, half-open (0-based) -> 1-start, fully-closed (1-based)，只需要在start那里+1，而end不用動
• 1-start, fully-closed (1-based) -> 0-start, half-open (0-based) ，只需要在start那里-1，而end不用動

一些轉(zhuǎn)換實例

USCS基因組瀏覽器的轉(zhuǎn)換

The UCSC Genome Browser Coordinate Counting Systems 文章中提到了，USCS內(nèi)源數(shù)據(jù)和基因組瀏覽器上數(shù)據(jù)的轉(zhuǎn)換

還提到了用liftOver軟件轉(zhuǎn)換的一些實例，大家可以自己去看

前面那篇biostar文章中，也提到了他們課題組自己開發(fā)的一個格式轉(zhuǎn)換工具convert_zero_one_based

IGV的格式轉(zhuǎn)換

親測narrowPeak放入IGV里面會start會加1，大家可以隨便拿個narrowPeak試試

Chr1    802633  802999 -> Chr1  802634  802999

Grange轉(zhuǎn)換成bed

可以參考Question: How To Write Data In A Granges Object To A Bed File. 這個問答

里面提到了手動轉(zhuǎn)換

# 注意 1 -> 0的話，start要-1

gr <- GRanges(seqnames = Rle(c("chr1", "chr2", "chr1", "chr3"), c(1, 3, 2, 4)),
  ranges = IRanges(1:10, end = 7:16, names = head(letters, 10)),
  strand = Rle(strand(c("-", "+", "*", "+", "-")), c(1, 2, 2, 3, 2)))

df <- data.frame(seqnames=seqnames(gr),
  starts=start(gr)-1,
  ends=end(gr),
  names=c(rep(".", length(gr))),
  scores=c(rep(".", length(gr))),
  strands=strand(gr))

write.table(df, file="foo.bed", quote=F, sep="\t", row.names=F, col.names=F)

也可以用 rtracklayer 包轉(zhuǎn)換。

narrowPeak用ChIPseeker注釋

其實這篇文章的一開始起意于Y叔這篇文章基因組位置讀進去竟然有一個位移，才讓我花了時間去好好探究下 0 和 1。經(jīng)過上面的理解，我們就會明白，之所以輸出的結(jié)果會比，輸入的bed文件或者narrowPeak文件，start部分要+1。就是因為bed是0-based，而Y叔的注釋Peak的主體是Grange格式，其是1-based格式。如果我們?nèi)ゼ毦縔叔的代碼的話，就會發(fā)現(xiàn)其在讀取Peak那個函數(shù)中，在start部分+1了。

image

但還有個小問題就是，由于Txdb對象的seqlevels有時候跟我們narrowPeak的seqlevels是不匹配的，比如一個是1，另一個是Chr1。所以我們時常會自己手動轉(zhuǎn)換GRange。但annotatePeak對于+1這一步實際上是發(fā)生在讀取Peak那一步的，所以你如果輸入的是自己準備的GRange，那么annotatePeak就不會+1了。這就意味著，如果你是自己用GRange函數(shù)手動轉(zhuǎn)換bed成的GRange或者由其他包產(chǎn)生的GRange格式，就要自己考慮下是否+1這個問題了。