大神一句話,，菜鳥跑半年,。我不是大神，但我可以縮短你走彎路的半年~
就像歌兒唱的那樣,，如果你不知道該往哪兒走,，就留在這學(xué)點(diǎn)生信好不好~
這里有豆豆和花花的學(xué)習(xí)歷程，從新手到進(jìn)階,，生信路上有你有我,！
豆豆寫于19.3.30
我們肯定都遇到許多利用坐標(biāo)去處理范圍信息的需求，比如要定位基因組的某個(gè)位置,，這個(gè)位置可能代表了gene model ,、genetic variants(包括了SNPs、inser‐tions/deletions),、transposable elements ,、binding sites；又或者想看看染色體某個(gè)區(qū)域的GC含量,、統(tǒng)計(jì)overlap、計(jì)算coverage,、提取序列等,。這些都屬于Range Data的處理范圍。本次著重看原理部分,。
推薦一本英文版的書 Bioinformatics Data Skills
處理Range Data一般有兩種途徑：R中的GenomicRanges和Linux中的bedtools
什么是RangesData,？ Ranges are integer intervals that represent a subsequence of consecutive positions on a sequence like a chromosome.
它是一種坐標(biāo)，記錄序列信息,?；蚪M的坐標(biāo)規(guī)定都是整數(shù)，所以不會(huì)出現(xiàn)類似這樣的坐標(biāo)(50,403,503.53)
坐標(biāo)有三大要素： 染色體/序列名稱 （因?yàn)橛械幕蚪M并沒(méi)拼接完,，還在scaffold或contig階段）：每個(gè)基因組都由一組染色體序列組成,，我們需要指定在哪個(gè)大范圍中。但是目前染色體命名沒(méi)有一個(gè)標(biāo)準(zhǔn),，比如UCSC和NCBI的染色體是chr開頭chr1,，而emsembl直接數(shù)字1
區(qū)間 ：比如114,414,997 to 114,693,772，由起始位點(diǎn)和終止位點(diǎn)組成
鏈 ：因?yàn)镈NA是雙鏈,，所以基因的特征信息(gene feature)可以存儲(chǔ)在正鏈(positive/forward)或者負(fù)鏈(negative/reverse)
需要注意的是,，坐標(biāo)信息與參考基因組相關(guān),，因此在討論ranges的時(shí)候要注意基因組的版本號(hào)，比如chr15:27,754,876-27,755,076這個(gè)區(qū)間在不同版本的基因組中就表示不同信息,，尤其是在和別人共享信息時(shí)可以說(shuō)我得到的這個(gè)范圍是基于GRCh38的,，或者基于GRCh37/hg19的
如果之前根據(jù)舊版本的基因組得到的坐標(biāo)要遷移到新版 ，可以使用一些工具,，比如CrossMap 【支持BED,、GFF/GTF、SAM/BAM,、Wiggle,、VCF文件格式在不同基因組版本之間的切換】、NCBI Genome Remapping Service 【是一個(gè)網(wǎng)頁(yè)工具】,、LiftOver 【基于UCSC基因組瀏覽器】
不同的range類型 下圖中x和y存在1個(gè)bp的overlap ,；z沒(méi)有任何overlap，只是自己跨越了1個(gè)bp,；
通過(guò)箭頭可以知道,，x和z都是正鏈，分別表示ACTT和C ,，y是在負(fù)鏈,，其堿基信息就要從3'向5'看，AGCCTTCG
兩套坐標(biāo)基準(zhǔn) 0-based：就是說(shuō)序列的第一位坐標(biāo)記為0,，最后一位坐標(biāo)是序列長(zhǎng)度-1,，遵循半開半閉區(qū)間[start,end) ，包含左邊不包含右邊,。例如：[1,5) 表示坐標(biāo)為1,，2，3,，4 ,，這個(gè)模式就和Python一樣
>>>'CTTACTTCGAAGGCTG'[1:5]
1-based：這個(gè)好像比較符合我們的習(xí)慣，從1開始,，遵循雙閉區(qū)間[start,end] ,，兩邊都包括。例如：[2,5] 就是2,3,4,5 ,，這個(gè)模式就和R一樣
> substr('CTTACTTCGAAGGCTG', 2, 5)這兩種系統(tǒng)各有優(yōu)劣：
盡管我們認(rèn)為1-based更符合自然計(jì)數(shù)法則,，但是有些情況下并不好用，比如：end-start 就好，這也比較好理解；但是1-based系統(tǒng)需要end-start+1
另外,，0-based系統(tǒng)支持zero-width feature ,，常用來(lái)描述兩個(gè)堿基之間的位置，比如在上圖中我們現(xiàn)在找一個(gè)酶切位點(diǎn)[12,12) ,，然后序列就被分成了CTTACTTCGAAGG 和CTG ,；這一點(diǎn)1-based系統(tǒng)也不能實(shí)現(xiàn)，它最小就是1個(gè)堿基
不同文件支持的坐標(biāo)系統(tǒng)不同：
比較麻煩的鏈信息 從https://www./p/3423/看了幾個(gè)重要的概念：
forward strand, this means reading left-to-right, and for the reverse strand it means right-to-left
A gene can live on a DNA strand in one of two orientations. The gene is said to have a coding strand (also known as its sense strand ), and a template strand (also known as its antisense strand ).
mRNA sequence always corresponds to the 5-3 coding sequence of a gene.
mRNA matches the coding sequence of the gene, not the template sequence(看圖https://en./wiki/File:Simple_transcription_elongation1.svg)
注意：上圖中的文件(除了Blast)都是基于參考基因組生成，而參考基因組序列是以正鏈為基準(zhǔn)
這也就是說(shuō)：我們從UCSC,、NCBI或Ensembl下載的參考基因組都是正鏈堿基序列,。
因此如果某個(gè)基因存在于負(fù)鏈 ,，從bed,、GTF等文件中看到的坐標(biāo)比如是chr1：2,473,087-2,492,258，抽取的這一區(qū)間序列比如是TCTTTAC...CCGAA ,，其實(shí)真正NCBI記錄的基因序列是TTCGG...GTAAAGA ,，與bed或GTF中記錄的序列正好是反向互補(bǔ)
因此，如果想從參考基因組中抽出位于負(fù)鏈的基因序列,，需要：1.先抽出參考基因組給的序列,；2.將序列反向互補(bǔ)。
對(duì)于轉(zhuǎn)錄本在負(fù)鏈的情況,，exon的實(shí)際位置也變了：原來(lái)在bed或GTF中forward 5'=>3'記錄的第一個(gè)位置,，實(shí)際上是在轉(zhuǎn)錄本的末尾；記錄的最后一個(gè)位置
所以,，我們?cè)贕TF文件中看到的正負(fù)鏈信息就十分有用了：
如果記錄+，表示在正鏈,，那么沒(méi)有問(wèn)題,，和文件中記錄的位置和序列都一樣；
如果記錄-,，那么真實(shí)信息一定是和文件記錄的信息反向互補(bǔ)的,，實(shí)際的位置也會(huì)改變
花了15分鐘思考做出來(lái)這個(gè)圖
最后的補(bǔ)充：
鏈的信息的確很重要，但是沒(méi)有做鏈特異性建庫(kù)時(shí)（表示我們是不知道鏈的方向信息的）,，如果要統(tǒng)計(jì)有多少reads比對(duì)到了某個(gè)特定的基因,，一般reads會(huì)在兩條鏈都有比對(duì)，那么這時(shí)一般會(huì)將兩條鏈的比對(duì)結(jié)果都統(tǒng)計(jì)上,，以免漏掉真正基因的區(qū)域
