測(cè)序數(shù)據(jù)比對(duì)和變異檢測(cè)

全基因組重測(cè)序是對(duì)已知基因組序列的物種進(jìn)行不同個(gè)體的基因組測(cè)序，并在此基礎(chǔ)上對(duì)個(gè)體或群體進(jìn)行差異性分析。所以首先我們需要某個(gè)物種的全基因組序列和該物種的某個(gè)個(gè)體的基因組序列。

重測(cè)序數(shù)據(jù)分析流程

重測(cè)序數(shù)據(jù)分析流程

一 、參考基因組的下載

進(jìn)入NCBI下載ASM1252v1的參考基因組

wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/012/525/GCA_000012525.1_ASM1252v1/GCA_000012525.1_ASM1252v1_genomic.fna.gz

gunzip GCA_000012525.1_ASM1252v1_genomic.fna.gz
下載地址的獲取

二、測(cè)試序列（即從測(cè)序數(shù)據(jù)）的準(zhǔn)備：

下載鏈接：
鏈接：https://pan.baidu.com/s/1cpLI4RqEs7_Ulkkl76w75g
提取碼：eg7g
復(fù)制這段內(nèi)容后打開百度網(wǎng)盤手機(jī)App，操作更方便哦

三、創(chuàng)建項(xiàng)目目錄

mkdir ~/Seqs/bwa_test

cp GCA_000012525.1_ASM1252v1_genomic.fna bwa_test/

cp test_* bwa_test/

bwa的介紹

  1. 即Burrows-Wheeler-Alignment Tool。
  2. 是一種能夠?qū)⒉町惗容^小的序列比對(duì)到一個(gè)較大的參考基因
    組上的軟件包。它由三個(gè)不同的算法:
    – BWA-backtrack: 是用來(lái)比對(duì) Illumina 的序列的,reads 長(zhǎng)度最長(zhǎng)
    能到 100bp。
    – BWA-SW: 用于比對(duì) long-read ,支持的長(zhǎng)度為 70bp-1Mbp;同時(shí)
    支持剪接性比對(duì)。
    – BWA-MEM: 和SW都支持較長(zhǎng)的read長(zhǎng)度,同時(shí)都支持剪接性比
    對(duì)(split alignments),但是BWA-MEM是更新的算法,也更快,
    更準(zhǔn)確,且 BWA-MEM 對(duì)于 70bp-100bp 的 Illumina 數(shù)據(jù)來(lái)說(shuō),
    效果也更好些。

bwa的安裝

sudo apt install bwa

bwa
bwa

bwa的使用

bwa的使用需要兩種輸入文件:
– 索引的Reference genome data(fasta格式 .fa, .fasta, .fna)
– Short reads data (fastaq格式 .fastaq, .fq)即重測(cè)序數(shù)據(jù)

四、參考基因組索引的建立

bwa index GCA_000012525.1_ASM1252v1_genomic.fna

五、reads比對(duì)到參考序列得到sam文件

bwa mem GCA_000012525.1_ASM1252v1_genomic.fna test_7942raw_1.fq.gz test_7942raw_2.fq.gz > test_bwa_7942.sam
2018-11-28 13-32-22屏幕截圖.png

SAM文件

1.SAM格式主要應(yīng)用于測(cè)序序列mapping到基因組上的結(jié)果表示,當(dāng)然也可以表示任意的多重比對(duì)結(jié)果。
2.SAM是一種序列比對(duì)格式標(biāo)準(zhǔn),由Heng Li (Sanger)制定,是以TAB為分割符的文本格式。
3.SAM的全稱是sequence alignment/map format。
sam格式詳解見https://en.wikipedia.org/wiki/SAM_%28file_format%29

less test_bwa_7942.sam
2018-11-28 13-37-46屏幕截圖.png

BAM格式

1.sam是帶有比對(duì)信息的序列文件(即告訴你這個(gè)reads在染色體上的位置等),用于儲(chǔ)存序列數(shù)據(jù)。
2.BAM 也是存儲(chǔ)序列比對(duì)信息的文件格式,但是是以二進(jìn)制存儲(chǔ),可以節(jié)約空間,計(jì)算機(jī)可讀。
3.二者都是fastq文件經(jīng)過(guò)序列比對(duì)或者mapping后輸出的格式;其儲(chǔ)存的信息都是一致的。

SAMtools

samtools是一個(gè)用于操作sam和bam文件的工具合集。

sudo apt install samtools

samtools

2018-11-30 11-19-02屏幕截圖.png

Official Document:http://www.htslib.org/doc/samtools.html

bam 與sam的格式轉(zhuǎn)換，bam格式以二進(jìn)制存儲(chǔ)文件，更加節(jié)約空間。

為參考基因組建立索引,生成了prefix.fai文件

samtools faidx GCA_000012525.1_ASM1252v1_genomic.fna

sam文件轉(zhuǎn)為bam文件

samtools view -bhS -t GCA_000012525.1_ASM1252v1_genomic.fna.fai -o test_bwa_7942.bam test_bwa_7942.sam

為bam文件排序,sort只能為bam文件排序,而不能為sam;不同版本samtools sort命令的-o參數(shù)不同

samtools sort test_bwa_7942.bam -o test_bwa_7942.bam.sorted

為排序的bam文件建立索引. *.bai文件

samtools index test_bwa_7942.bam.sorted

samtools tview:顯示reads比對(duì)基因組的情況

samtools tview test_bwa_7942.bam.sorted GCA_000012525.1_ASM1252v1_genomic.fna
2018-11-30 11-31-34屏幕截圖.png

測(cè)試參考基因組每個(gè)位點(diǎn)或一段區(qū)域的測(cè)序深度(不是常說(shuō)的平均測(cè)序深度)

samtools depth test_bwa_7942.bam.sorted >depth.txt
2018-11-30 11-34-07屏幕截圖.png

samtools flagstat:統(tǒng)計(jì)比對(duì)結(jié)果

samtools flagstat test_bwa_7942.bam.sorted
2018-11-30 11-35-37屏幕截圖.png

-總共的reads數(shù)(QC-passed or failed)
-總體上reads的匹配率
-有多少reads是屬于paired reads
-reads1中的reads數(shù)
-reads2中的reads數(shù)
-properly paired:正確配對(duì)的reads數(shù)量
-with itself and mate mapped:一對(duì)reads均比對(duì)上的reads數(shù)
-singletons:只有單條reads比對(duì)上的reads數(shù)

samtools rmdup:去除PCR重復(fù)

samtools rmdup test_bwa_7942.bam.sorted output.bam

samtools mpileup:生成bcf文件

該命令用于生成bcf文件,再使用bcftools進(jìn)行SNP和Indel的分析

samtools mpileup -gf GCA_000012525.1_ASM1252v1_genomic.fna test_bwa_7942.bam.sorted > test_bwa_7942.bcf

mpileup生成的結(jié)果是pileup格式包含6行:

  1. 參考序列名;
  2. 位置;
  3. 參考?jí)A基;
  4. 比對(duì)上的reads數(shù);
  5. 比對(duì)情況;
    – ‘.’代表與參考序列正鏈匹配。
    – ‘,’代表與參考序列負(fù)鏈匹配。
    – ‘ATCGN’代表在正鏈上的不匹配。
    – ‘a(chǎn)tcgn’代表在負(fù)鏈上的不匹配。
    – ‘*’代表模糊堿基
    – ‘’代表匹配的堿基是一個(gè)read的開始;’’后面緊跟的ascii碼減去33代表比對(duì)質(zhì)量;這兩個(gè)符號(hào)修飾的是后面的堿基,其后緊跟的堿基(.,ATCGatcgNn)代表該read的第一個(gè)堿基。
    – ‘$’代表一個(gè)read的結(jié)束,該符號(hào)修飾的是其前面的堿基。
    – 正則式’+[0-9]+[ACGTNacgtn]+’代表在該位點(diǎn)后插入的堿基;比如上例中在scaffold_1的2847后插入了9個(gè)長(zhǎng)度的堿基acggtgaag。表明此處極可能是indel。
    –正則式’-[0-9]+[ACGTNacgtn]+’代表在該位點(diǎn)后缺失的堿基;
    6.比對(duì)上的堿基的質(zhì)量

六、基因變異檢測(cè)軟件

從bcf文件(samtools mpileup)中檢測(cè)基因變異位點(diǎn): bcftools，從bam文件中檢測(cè)基因變異位點(diǎn):GATK，基因變異的主要類型:SNP, indel等。bcftools是samtools中附帶的軟件,在samtools的安裝文件夾中可以找到。

bcftools call -vm test_bwa_7942.bcf -o test_bwa_7942.variants.bcf

bcftools view -v snps,indels test_bwa_7942.variants.bcf>test_bwa_7942.snps.vcf

less test_bwa_7942.snps.vcf

變異位點(diǎn)的過(guò)濾——bcftools filter

bcftools filter -o test_bwa_7942.snps.filtered.vcf -i 'QUAL>20 &&DP>5' test_bwa_7942.snps.vcf

這里簡(jiǎn)單過(guò)濾:QUAL小于等于20,DP值小于等于5的位點(diǎn),-i – include 只保留滿足該條件的位點(diǎn)。

七、變異基因注釋軟件

Annovar安裝和運(yùn)行

wget http://www.openbioinformatics.org/annovar/download/0wgxR2rIVP/annovar.latest.tar.gz

tar -zvxf annovar.latest.tar.gz -C ~/BioSofts/

echo 'export PATH=~/BioSofts/annovar:$PATH'>>~/.bashrc

source ~/.bashrc

Annovar三種注釋模式

  1. gene-based annotation:判斷SNV或CNV是否造成蛋白編碼或氨基酸的改變,可用基因命名系統(tǒng)包括RefSeq,UCSC,ENSEMBL,GENCODE, AceView等。
  2. region-based annotation:變異位于染色體哪個(gè)區(qū)域,預(yù)測(cè)轉(zhuǎn)錄因子結(jié)合位點(diǎn)、SD區(qū)域等
  3. filter-based annotation:鑒定在特定數(shù)據(jù)庫(kù)中記錄的變異,如是否在dbSNP中被報(bào)道。
    http://annovar.openbioinformatics.org/en/latest/

Annovar主要程序和目錄

-annotate_variation.pl #主程序,功能包括下載數(shù)據(jù)庫(kù),三種不同的注釋
-coding_change.pl #可用來(lái)推斷蛋白質(zhì)序列
-convert2annovar.pl #將多種格式轉(zhuǎn)為.avinput的程序
-retrieve_seq_from_fasta.pl #用于自行建立其他物種的轉(zhuǎn)錄本
-table_annovar.pl #注釋程序,可一次性完成三種類型的注釋
-variants_reduction.pl #可用來(lái)更靈活地定制過(guò)濾注釋流程
-example/ #存放示例文件
-humandb/ #人類注釋數(shù)據(jù)庫(kù)

生成annovar輸入文件

convert2annovar.pl -format vcf4 test_bwa_7942.snps.filtered.vcf >test_bwa_7942.snps.avinput

Avinput文件,每行代表一個(gè)位點(diǎn),前5列依次為:
– chromosome
– start position
– end position
– the reference nucleotides
– the observed nucleotides
– 其余列:可有可無(wú);如果有,在輸出文件中會(huì)原樣輸出。

注釋變異基因位點(diǎn)

annotate_variation.pl --geneanno --dbtype refGene --buildver 7942-genome-new test_bwa_7942.snps.avinput ~/BioSofts/annovar/humandb/

--geneanno: 采用gene-based annotation注釋模式;
--dbtype :數(shù)據(jù)庫(kù)為refGene;還有knowGene,ensGene等
--buildver:為參考基因組版本
最后為數(shù)據(jù)庫(kù)所在目錄
生成avinput.variant_function和avinput.exonic_variant_function后綴的兩個(gè)結(jié)果文件。

Annovar結(jié)果解析

1.avinput.variant_function文件
2.包括所有突變信息:
3.第一列:variant effects, 變異分類,如intergenic, intronic, non-synonymous SNP, frameshift deletion, large-scale duplication等
4.第二列:基因名或Symbol
5.其余列:與avinput輸入文件相同:染色體、start、end、ref_nt、obs_nt等

需要下載或自定義注釋數(shù)據(jù)庫(kù)

1.官方提供了一些基因組注釋數(shù)據(jù)庫(kù)
http://annovar.openbioinformatics.org/en/latest/user-guide/download/
2.如果沒(méi)有,這需要自定義注釋數(shù)據(jù)庫(kù)

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容