宅在家兩個多月戒祠，不知不覺已經是春天了骇两，也許距離返校的日子更近了吧...

自己拍的

變異和突變

變異，指的是實際測序數據與國際規(guī)定的參考基因組之間的區(qū)別姜盈。很多變異其實只是造成人類多樣性的原因低千。突變，指的是那些與疾病相關的變異馏颂。
舉個例子：ENSEMBL等規(guī)定的人類參考基因組文件某位置是AAAAA示血，然后一個人實際測序得到的序列為AGCAA，那么相比于參考基因組救拉，這個人就有2個變異位點矾芙。對于第2個位置，如果查看所有已知的測序近上，絕大部分人都是G，說明是參考基因組出現了問題拂铡，這個變異就不能稱作突變壹无。對于第3個位置，如果查看所有已知的測序感帅，絕大部分人都是A斗锭，而恰好有一個人不是A，但他是個患者失球，那么這個變異就是突變了岖是。

變異的類型

SNP（single nucleotide polymorphism）：單核苷酸多態(tài)性。個體間基因組DNA序列同一位置單個核苷酸變異（替換实苞、插入或缺失）所引起的多態(tài)性豺撑。在人類基因組中SNP分布普遍并且密度較大，總數超過107黔牵， 平均每300bp（也有說1kbp）就有一個SNP聪轿。或稱單核苷酸位點變異SNV猾浦。
INDEL（insertion-deletion）：插入和缺失陆错。基因組上小片段（>50bp）的插入或缺失。
CNV（copy number variation）：基因組拷貝數變異金赦。基因組中大片段的DNA形成非正常的拷貝數量音瓷。比如一個基因在染色體的一條染色單體上的數目為1，但是在染色體復制過程中夹抗，復制結束后該基因在染色單體數目由1變成了2或者n绳慎。它發(fā)生的頻率遠遠高于染色體結構變異，并且整個基因組中覆蓋的核苷酸總數大大超過SNP的總數。
SV（structure variation）：結構變異偷线。染色體大片段的插入與缺失磨确，染色體內部的某區(qū)域發(fā)生翻轉顛換，兩條染色體之間發(fā)生重組声邦。

圖片來源于網絡

SNP

一般情況下只分析SNP乏奥，其它類型的變異分析有難度或不準確。
來自兩個不同個體的DNA片段AAGCCTA和AAGCTTA為等位基因亥曹。幾乎所有常見的SNP位點只有兩個等位基因邓了。
在人體中，SNP的發(fā)生機率大約是0.1%媳瞪，也就是每1000個堿基對就可能有一個SNP（密度高）骗炉。對疾病發(fā)生和藥物治療有重大影響的SNP，估計只占數以百萬計SNP的很小一部分蛇受。
SNP位點的分布是不均勻的句葵，在非轉錄序列比在轉錄序列更常見。編碼區(qū)的單核苷酸多態(tài)性——編碼 SNP（coding SNP兢仰，cSNP）也有同義和非同義兩種類型乍丈，非同義SNP會改變蛋白質的氨基酸序列“呀基因非編碼區(qū)轻专、基因間隔區(qū)的SNP仍然可能影響轉錄因子結合、剪接等過程察蹲。
從演化的觀點來看请垛，SNP具有相當程度的穩(wěn)定性，即使經過代代相傳洽议，SNP所引起的改變卻不大宗收，因此可用以研究族群演化。

基于測序數據檢測SNP的方法

圖片來源于網絡

基于測序數據檢測結構變異SV的方法

圖片來源于網絡

以HISAT2+Samtools+bcftools為例亚兄，識別變異的流程

HISAT2

HISAT2 是一款利用改進的BWT算法進行序列比對的軟件镜雨。由約翰霍普金斯大學計算生物學中心（CCB at JHU）開發(fā)，是TopHat的升級版本儿捧，速度提高了50倍荚坞。利用 HISAT2 + StringTie 流程，可以快速地分析轉錄組測序數據菲盾，獲得每個基因和轉錄本的表達量颓影。

首先需要構建參考基因組索引用于下一步的比對。HISAT2提供了兩個腳本用于從基因組注釋GTF文件中提取剪接位點和外顯子位置懒鉴，基于這些特征诡挂，可以使 RNA-Seq reads 比對更加準確碎浇。然后再進行reads mapping。

#（1）下載參考基因組和注釋文件
# 例如：dmel.genome.fa.gz和dmel.annotation.gtf.gz
wget  dmel.genome.fa.gz
gzip  -d  dmel.genome.fa.gz  #或gunzip  dmel.genome.fa.gz
#（2）提取外顯子和剪接位點
extract_splice_sites.py  dmel.annotation.gtf  >  dmel.ss
extract_exons.py  dmel.annotation.gtf  >  dmel.exon
extract_snps.py  snp142Common.txt  >  genome.snp
#（3）構建參考基因組索引
hisat2-build  --ss  dmel.ss  --exon  dmel.exon  (--snp  genome.snp)  dmel.genome.fa  dmel
# 得到8個索引文件璃俗，以dmel為共同的前綴奴璃，dmel.1~8.ht2
#（4）也可以下載現成的index
wget  bdgp6.tar.gz
tar  -zxvf  bdgp6.tar.gz
#（5）reads mapping
hisat2  -p  8  (--dta)  -x  dmel  -1  sample_1.fq.gz  -2  sample_2.fq.gz  -S  sample.sam  (&)
# hisat2  [options]*  -x <ht2-idx>  {-1 <m1>  -2 <m2> | -U <r> | --sra-acc <SRA accession number>}  [-S <sam>]
# -p  線程數
# -t  記錄時間
# --dta   輸出專門為轉錄本組裝的比對結果（如果比對結果后續(xù)需要使用StringTie進行轉錄本組裝，則需要加入--dta選項）
# -x  指定基因組索引的文件前綴
# -S  指定輸出的SAM文件城豁。默認輸出到標準輸出苟穆，比對結束后統(tǒng)計結果輸出到標準錯誤輸出。
# -f  表示輸入文件格式為fasta（默認）唱星，-q表示輸入文件格式為fastq雳旅。可以是gzip壓縮的文件间聊。
# -phred33  輸入的FASTQ文件堿基質量值編碼標準為phred33（默認）攒盈。

比對結果：

比對輸出內容.PNG

Samtools

SAM（sequence Alignment/mapping）數據格式是目前高通量測序中存放比對數據的標準格式。BAM是SAM的二進制格式哎榴。使用samtools將sam文件轉化為bam文件型豁，并進行排序。

#（1）將sam文件轉化為bam文件
samtools view  -bS  sample.sam  >  sample_unsorted.bam
# -b：輸出為bam格式
# -S：出入為sam格式
#（2）將bam文件進行排序
samtools sort  -@  8  sample_unsorted.bam  -o  sample.sorted.bam
# -o：輸出文件的名字
# -@：線程數
#默認按照染色體位置進行排序尚蝌，-n是根據read名進行排序偷遗，-t 根據TAG進行排序。
#（3）版本1.4后可直接簡化
samtools sort  -@  8  -o  sample.sorted.bam  sample.sam
#（4）索引
samtools index  sample.sorted.bam
#得到sample.sorted.bam.bai索引文件

SAM文件：

SAM文件.PNG

bcftools

vcf格式（Variant Call Format）是存儲變異位點的標準格式驼壶，用于記錄variants（SNP / InDel）。BCF是VCF的二進制文件喉酌。

#（1）SNP calling
# mpileup命令：得到染色體上每個堿基的比對情況的匯總
bcftools mpileup  -Ob  -o  sample.bcf  -f  dmel.genome.fa  sample.sorted.bam
# 輸入BAM文件sorted.bam
# -f  --fasta-ref  指定參考序列的fasta文件
# -O  指定輸出文件的類型热凹，壓縮的BCF(b)，未壓縮的BCF(u)泪电，壓縮的VCF(z)般妙，未壓縮的VCF(v) 
# -o  指定輸出文件的名字sample.bcf
# call命令：執(zhí)行SNP calling 
bcftools call  -vmO  z  -o  sample.vcf.gz  sample.bcf
# -v  只輸出變異位點的信息，如果一個位點不是snp/indel則不會輸出
# 兩種calling算法：-c參數對應consensus-caller算法相速，-m參數對應multiallelic-caller算法碟渺，后者更適合多種allel和罕見變異的calling
#（2）tabix
tabix  -p  vcf  sample.vcf.gz
# 輸入為壓縮文件vcf.gz，生成的索引文件為sample.vcf.gz.tbi
#（3）index對vcf文件建立索引
bgzip  view.vcf  #輸入的VCF文件必須是bgzip壓縮后的文件
gunzip  view.vcf.gz  #解壓縮
bcftools index  view.vcf.gz  #生成索引文件view.vcf.gz.csi
bcftools index  -t  view.vcf.gz  #生成索引文件view.vcf.gz.tbi
#（4）query通過表達式來指定輸出格式
bcftools query  -f  '%CHROM\t%POS\t%REF\t%ALT[\t%SAMPLE=%GT]\n'  view.vcf.gz
# -f  通過一個表達式來指定輸出格式
# %CHROM  代表VCF文件中染色體那一列突诬，其他的列POS, ID, REF, ALT, QUAL, FILTER也是類似的寫法
# []  對于FORMAT字段的信息用中括號括起來
# %SAMPLE  代表樣本名稱
# %GT  代表FORMAT字段中genotype的信息
# \t  代表制表符分隔
# \n  代表新的一行
#（5）sort按照染色體位置進行排序
bcftools sort  view.vcf.gz  -o  sort.view.vcf
#（6）filter過濾不可靠位點
bcftools filter  -O  z  -o  sample_filtered.vcf.gz  -s  LOWQUAL -i'%QUAL>10'  sample.vcf.gz
# -O  --output-type 輸出的格式苫拍，z和v都行，壓縮的VCF(z)旺隙，未壓縮的VCF(v)
# -o  --output 輸出文件的名稱
# -s  --soft-filter 將過濾掉的位點用字符串注釋
#（7）view命令用于VCF和BCF格式的轉換
bcftools view  view.vcf.gz  -O  u  -o  view.bcf  
bcftools view  view.vcf.gz  -s  NA00001,NA00002  -o  subset.vcf 
bcftools view  view.vcf.gz  -k  -o  known.vcf 
# -O  指定輸出文件的類型绒极，壓縮的BCF(b)，未壓縮的BCF(u)蔬捷，壓縮的VCF(z)垄提，未壓縮的VCF(v) 
# -o  指定輸出文件的名字
# -s  想要保留的樣本信息榔袋，多個樣本用逗號分隔；如果樣本名稱添加了^前綴铡俐，代表去除這些樣本凰兑，比如-s ^NA00001,NA00002
# -k  表示篩選已知的突變位點，即ID那一列值不是.的突變位點
#（8）stats命令用于統(tǒng)計VCF文件的基本信息
bcftools stats  view.vcf  >  view.stats
bcftools stats  -F  dmel.genome.fa  -s  -  sample.vcf.gz  >  sample.vcf.gz.stats 
# -F  --fasta-ref  faidx indexed reference sequence（參考基因組） file to determine INDEL context
# -s  list of samples for sample stats, "-" to include all samples 統(tǒng)計的樣本列表审丘，-代表所有樣本
#（9）plot-vcfstats命令進行可視化
plot-vcfstats  sample.vcf.gz.stats  -p  plots/sample.vcf.gz.stats
# -p  指定輸出結果的目錄
#（這個腳本位于bcftools安裝目錄的misc目錄下吏够，依賴latex 生成pdf 文件。這里需要matplotlib）

stats統(tǒng)計文件：

stats結果.PNG