劉小澤寫于18.8.10
老規(guī)矩谨读，先上官網(wǎng)https://samtools.github.io/bcftools/bcftools.html
也許你會(huì)問局装，call是什么鬼？不能通俗易懂點(diǎn)嗎劳殖？好的～我的理解是召喚铐尚，BCFtools就是你手中的魔法棒，揮一揮哆姻，召喚變異位點(diǎn)
學(xué)一項(xiàng)工具宣增，掌握一種技術(shù)的目的是為了能夠更好地發(fā)現(xiàn)并解決問題，而不是為了跑流程而跑流程矛缨，所有的數(shù)據(jù)分析都要建立在個(gè)人背景知識(shí)儲(chǔ)備和對(duì)實(shí)際項(xiàng)目的認(rèn)知上爹脾。因此今天，我也只是用這一個(gè)工具來做一個(gè)引子箕昭，從變異開始學(xué)起灵妨。
要學(xué)習(xí)任何事物，都要從“是什么落竹？為什么泌霍？怎么做？“入手述召，那么～
變異是什么
你我之間朱转，只差1%的不同，而這點(diǎn)不同桨武，卻造就了幾十億人口的千差萬別
變異是相對(duì)的肋拔，是在彼此的比較中產(chǎn)生的。我們檢測(cè)基因組變異呀酸，都需要參考基因組，對(duì)于人類而言琼梆，一般就是用hg19性誉、hg38窿吩、b37版本。但是參考基因組也是選志愿者測(cè)的错览，并不完美纫雁，因此也不能百分之百就說參考基因組變異就是0，所有的都要參考這個(gè)基因組來確定倾哺。
如果只有這一個(gè)條件轧邪，不足以讓我們信服。因?yàn)槿绻麢z測(cè)出來一下假陽性的變異位點(diǎn)羞海，后果可能病急亂投醫(yī)忌愚。因此，在做全基因組檢測(cè)過程中却邓，還需要選用其他幾個(gè)大樣本群體的變異位點(diǎn)數(shù)據(jù)庫進(jìn)行對(duì)照硕糊，盡可能排除那些假陽性變異。
基因組變異一般包括：
- 單堿基變異腊徙，學(xué)名單核苷酸多態(tài)性（SNP）简十，原來的定義是單個(gè)堿基導(dǎo)致的群體中廣泛存在的（約1%）的多態(tài)性，后來指與參考基因組不同的位點(diǎn)撬腾，它最常見也最簡單螟蝙；
- 小片段的插入與缺失（合稱InDel），一般發(fā)生在基因組上短的有序的基因片段民傻，長度小于50bp
- 更大范圍的結(jié)構(gòu)性變異（SV）胰默，研究的熱門，長度大于50bp的片段的插入饰潜、缺失（Big Indel）初坠，染色體倒位（Inversion）、染色體之間或者內(nèi)部發(fā)生易位（Translocation）彭雾、拷貝數(shù)變異（CNV）碟刺、串聯(lián)重復(fù)（Tandem repeat）、嵌合體（chimera）
為什么要找變異
人類基因組的變異與人類起源發(fā)展薯酝、疾病防控等密切相關(guān)半沽，二代測(cè)序的低成本帶來了很大的便利，讓我們有能力去做這件事吴菠；
-
變異不等于突變者填，有的變異缺失能夠引起疾病，比如75%的癌癥患者都存在TP53基因突變;TP53變異與近一半以上的癌癥發(fā)生有關(guān)做葵，包括肺癌占哟、胃癌、肝癌、膀胱癌榨乎、食管癌怎燥、乳腺癌、宮頸癌等多種癌癥蜜暑。但是還有一些基因铐姚，比如球場(chǎng)彎道超車的球員就是11號(hào)染色體上的ACTN3基因上MA000028位點(diǎn)發(fā)生變異。找變異肛捍，可以解釋更多生物現(xiàn)象
img
- 結(jié)構(gòu)變異的尋找更為重要隐绵，因?yàn)樗袝r(shí)會(huì)與環(huán)境互作增加出生缺陷、癌癥的風(fēng)險(xiǎn)
如何尋找變異
- 將reads比對(duì)到參考基因組
- 矯正比對(duì)（可選）
- 從比對(duì)結(jié)果中確定變異位點(diǎn)
- 根據(jù)某些需求過濾
- 變異位點(diǎn)注釋
最常用的“魔法棒”
- bcftoolshttp://www.htslib.org/doc/bcftools.html
- FreeBayes: https://github.com/ekg/freebayes
- GATK: https://software.broadinstitute.org/gatk/
- VarScan2: http://varscan.sourceforge.net/
- 當(dāng)然還有根據(jù)體細(xì)胞拙毫、生殖細(xì)胞依许、家系數(shù)據(jù)優(yōu)化的軟件，比如GATK就針對(duì)體細(xì)胞恬偷、生殖細(xì)胞做了不同的流程悍手，目前算得上是找變異的黃金搭檔
根據(jù)這一篇跑一個(gè)小流程，加深理解袍患。唉坦康，現(xiàn)在這學(xué)習(xí)，不跑個(gè)數(shù)據(jù)诡延，就像沒學(xué)習(xí)一樣
第一步 下載數(shù)據(jù)
需要用到EMBOSS工具包滞欠，詳細(xì)介紹https://wenku.baidu.com/view/5bdf981c55270722192ef7e3.html?pn=NaN
官網(wǎng)：https://www.ebi.ac.uk/seqdb/confluence/display/THD/EMBOSS+seqret
#軟件安裝Emboss-6.6.0
cd ~/.biosoft
wget ftp://ftp.ebi.ac.uk/pub/software/unix/EMBOSS/EMBOSS-6.6.0.tar.gz
#下載參考數(shù)據(jù)
mkdir -p ~/project/ebola/ref
cd ~/project/ebola/ref
ACC=AF086833
REF=ref/$ACC.fa
efetch -db=nuccore -format=fasta -id=$ACC > $REF
bwa index $REF
#下載測(cè)試數(shù)據(jù)
mkdir -p ~/project/ebola/raw
cd ~/project/ebola/raw
fastq-dump -X 100000 --split-files SRR1553500 #取前10萬條
# 27M Aug 10 22:58 SRR1553500_1.fastq
# 27M Aug 10 22:58 SRR1553500_2.fastq
服務(wù)器100M獨(dú)立光纖還不如自己電腦下載速度快【兩個(gè)同時(shí)下載的】
第二步 比對(duì)
#學(xué)習(xí)一下教程中的shell腳本寫法
ACC=AF086833
REF=~/project/ebola/ref/$ACC.fa
SRR=SRR1553500
BAM=$SRR.bam
R1=${SRR}_1.fastq
R2=${SRR}_2.fastq
TAG="@RG\tID:$SRR\tSM:$SRR\tLB:$SRR"
bwa mem -R $TAG $REF $R1 $R2 | samtools sort > $BAM
samtools index $BAM
# 12M Aug 10 2018 SRR1553500.bam
# 200 Aug 10 23:04 SRR1553500.bam.bai
第三步 召喚變異
ACC=AF086833
REF=~/project/ebola/ref/$ACC.fa
samtools mpileup -bvf $REF SRR1553500.bam | \
bcftools call --ploidy 1 -vm -Ov > bcftools.vcf
未完待續(xù)。肆良。筛璧。
