重測序分析

重測序：是對已對已知基因組的物種進行測序，去挖掘不同個體和群體之間的差異性。

重測序分析內(nèi)容：

SNP，INDEL, SV , SNV

進化分析，群體結(jié)構(gòu)：

LD,FST


分析數(shù)據(jù)和流程。

測序數(shù)據(jù)：fastq 格式文件

序列比對軟件：BWA 軟件（快速把小片段比對到基因組上）


分析過程：

1.測序數(shù)據(jù)質(zhì)控(fastqc)

-t? 測序數(shù)據(jù)線程數(shù)，數(shù)目越多速度越快

-o 指定輸出目錄

$/share/nas2/genome/biosoft/fastqc/v0.10.1/fastqc -o ./fastqcout/ ./test-data/reads/test_1.fq ./test-data/reads/test_2.fq


2.序列比對（bwa）

三種比對方法,MAMA算法支持100bp以上比對

建立索引，幫助序列比對快速閱讀

bwa -index

/share/nas2/genome/biosoft/bwa/0.7.17/bwa/bwa index ../test-data/genome/cucumber_ChineseLong_v2_genome.fa

實際進行比對

/share/nas2/genome/biosoft/bwa/0.7.17/bwa/bwa mem ../test-data/genome/cucumber_ChineseLong_v2_genome.fa ../test-data/reads/test_1.fq ../test-data/reads/test_2.fq -o ./test.sam?

輸出sam文件

將SAM文件轉(zhuǎn)換為BAM文件（節(jié)省空間）

/share/nas2/genome/biosoft/samtools/current/samtools view ./bwa/test.sam -o ./bwa/test.bam

（查看比對bam文件samtools view ./bwa/test.bam |less -s）

(查看bam文件比對結(jié)果/share/nas2/genome/biosoft/samtools/current/samtools view ./bwa/test.sam -o ./bwa/test.bam)


將bam文件按照比對到參考基因組的順序進行排序：

/share/nas2/genome/biosoft/samtools/current/samtools sort ./bwa/test.bam -o ./bwa/test.sort.bam

(查看測序深度/share/nas2/genome/biosoft/samtools/current/samtools depth ./bwa/test.sort.bam)


/share/nas2/genome/biosoft/samtools/current/samtools faidx ./test data/genome/cucumber_ChineseLong_v2_genome.fa? (samtools 建立索引為后續(xù)突變位點做準(zhǔn)備，生成.fai文件)

/share/nas2/genome/biosoft/samtools/1.9/bin/samtools dict -o ./test-data/genome/cucumber_ChineseLong_v2_genome.dict ./test-data/genome/cucumber_ChineseLong_v2_genome.fa(samtools dict 生成基因組字典文件用于后續(xù)比對，生成.dict文件）

/share/nas2/genome/biosoft/samtools/1.9/bin/samtools index ./bwa/test.sort.bam(針對比對結(jié)果生成索引文件，生成。bai格式文件)


SNP calling：

1.標(biāo)記重復(fù)序列（由于PCR擴增或者錯誤測序，進行標(biāo)記）

picard和GATK都可以

picard 的MarkDuplicates可以標(biāo)記重復(fù)序列（用法，輸入排序好的bam文件，輸出標(biāo)記重復(fù)的bam文件，-M矩陣文件(后期需要讀取)，系統(tǒng)設(shè)置?MAX_OPTICAL_DUPLICATE_SET_SIZE 最大可選大小，一般512）

java -jar /share/nas2/genome/biosoft/picard/picard.jar MarkDuplicates MAX_OPTICAL_DUPLICATE_SET_SIZE=512 I= ./bwa/test.sort.bam O= ./duplication

_marking/dedup.bam M= ./duplication_marking/dedup.metrics

對生成的bam建立索引

/share/nas2/genome/biosoft/samtools/current/samtools index ./duplication_marking/dedup.bam


2.統(tǒng)計中間文件

java -jar /share/nas2/genome/biosoft/GATK/3.3-0/GenomeAnalysisTK.jar -T RealignerTargetCreator -R ./test-data/genome/cucumber_ChineseLong_v2_genome.fa -I ./duplication_marking/test_dedup.bam -o ./duplication_marking/test_dedup.bam_realign_interval

3.發(fā)生錯配的地方進行重新比對提高假陽性和假陰性。

java -jar /share/nas2/genome/biosoft/GATK/3.3-0/GenomeAnalysisTK.jar -T IndelRealigner -R ./test-data/genome/cucumber_ChineseLong_v2_genome.fa -targetIntervals ./duplication_marking/test_dedup.bam_realign.intervals -o ./localalign/test_relign.dedup.bam -I ./duplication_marking/test_dedup.bam(targetIntervals標(biāo)記重復(fù)序列不需要比對，對文件名敏感，必須是.intervals)


SNP位點

生成gvcf文件

--indelSizeToEliminateInRefModel 指定最大Indel 長度

--emitRefConfidence GVCF 指定生成文件的類型

--variant_index_type LINEAR 指定使用的索引模型

--variant_index_parameter 指定使用的索引策略

--sample_ploidy 指定分析物種的倍性

-nct 指定分析使用的CPU 核心數(shù)量

-o 指定輸出文件

-L 指定染色體位置文件

-I 指定輸入BAM 文件

java -jar /share/nas2/genome/biosoft/GATK/3.3-0/GenomeAnalysisTK.jar -T HaplotypeCaller -R ./test-data/genome/cucumber_ChineseLong_v2_genome.fa --indelSizeToEliminateInRefModel 50 --emitRefConfidence GVCF --variant_index_type LINEAR --variant_index_parameter 128000 --sample_ploidy 2 -nct 4 -o ./variant_calling/test_g.vcf? -I ./localalign/test_relign.dedup.bam


###合并gcvf文件

-T 指定使用的功能大類

-R 指定參考基因組

--disable_auto_index_creation_and_locking_when_reading_rods 在讀取時禁用自動索引創(chuàng)建和鎖定(防止出錯)

-o 指定輸出文件

--variant 指定輸入GVCF 文件

java -jar /share/nas2/genome/biosoft/GATK/3.3-0/GenomeAnalysisTK.jar -T CombineGVCFs -R ./test-data/genome/cucumber_ChineseLong_v2_genome.fa --disable_auto_index_creation_and_locking_when_reading_rods -o ./variant_calling/cohort.test.g.vcf --variant ./variant_calling/test_g.vcf

##將GVCF文件轉(zhuǎn)化為vcf文件

-T 指定使用的功能大類

-nt 指定使用的CPU 核心數(shù)量

-R 指定參考基因組

--disable_auto_index_creation_and_locking_when_reading_rods 在讀取時禁用自動索引創(chuàng)建和鎖定

-o 指定輸出文件

--variant 指定輸入文件

java -jar /share/nas2/genome/biosoft/GATK/3.3-0/GenomeAnalysisTK.jar -T GenotypeGVCFs -R ./test-data/genome/cucumber_ChineseLong_v2_genome.fa --disable_auto_index_creation_and_locking_when_reading_rods -o ./variant_calling/cohort.test.snp.indel.vcf --variant ./variant_calling/cohort.test.g.vcf

###合并VCF文件

java -cp /share/nas2/genome/biosoft/GATK/3.3-0/GenomeAnalysisTK.jar org.broadinstitute.gatk.tools.CatVariants -R ./test-data/genome/cucumber_ChineseLong_v2_genome.fa -out ./variant_calling/raw.vcf -assumeSorted -V ./variant_calling/cohort.test.snp.indel.vcf? ? ? ? ?


SNP 過濾

1.一般來說如果臨近的堿基出現(xiàn)SNP一般說來是存在的，一般來說5個堿基以內(nèi),用10個bp進行劃窗分析

perl /share/nas2/genome/biosoft/bcftools/bin/vcfutils.pl varFilter -w 5 -W 10 ./variant_calling/raw.vcf >./variant_calling/raw.vcf.tmp? ? ? ? ? ? ? ? ? ? ? ? ? ??

2.按照質(zhì)量值、位點深度、Fisher 檢驗值、比對質(zhì)量對SNP 位點進行過濾

-T 指定使用的功能大類

-R 指定參考基因組

-V 指定輸入文件

--filterExpression "QUAL<30||QD < 2.0 || FS > 60.0 || MQ < 40.0"? 指定過濾條件

--clusterWindowSize 指定聚類窗口大小

--clusterSize 指定聚類數(shù)量

--filterName my_snp_filter 指定過濾條件名體現(xiàn)在輸出文件中

-o raw.filter.vcf.tmp 指定輸出文件

java -jar /share/nas2/genome/biosoft/GATK/3.3-0/GenomeAnalysisTK.jar -T VariantFiltration -R ./test-data/genome/cucumber_ChineseLong_v2_genome.fa -V ./variant_calling/raw.vcf.tmp --filterExpression "QUAL<30||QD < 2.0 || FS > 60.0 || MQ < 40.0" --clusterWindowSize 5 --clusterSize 2 --filterName my_snp_filter -o ./variant_calling/raw.filter.vcf.tmp

3.保留通過的列

awk '$7 =="PASS"|| $0~/#/{print $0}' ./variant_calling/raw.filter.vcf.tmp > ./variant_calling/raw.filter.vcf


###拆分Indel和SNP文件

僅僅保留INDEL位點的文件

java -jar /share/nas2/genome/biosoft/GATK/3.3-0/GenomeAnalysisTK.jar -T SelectVariants -R ./test-data/genome/cucumber_ChineseLong_v2_genome.fa -V ./variant_calling/raw.filter.vcf -selectType INDEL -o ./variant_calling/raw.filter.indel.vcf? ? ? ?

僅僅保留SNP位點的文件

java -jar /share/nas2/genome/biosoft/GATK/3.3-0/GenomeAnalysisTK.jar -T SelectVariants -R ./test-data/genome/cucumber_ChineseLong_v2_genome.fa -V ./variant_calling/raw.filter.vcf -selectType SNP -o ./variant_calling/raw.filter.snp.vcf


###使用SNPEFF 對SNP 位點進行注釋

1.構(gòu)建SNPEFF 物種數(shù)據(jù)庫

data 的目錄不能更改，必要要有，下一級目錄取物種名字，添加兩個文件，基因組的gff和fa文件，必須命名為sequences.fa 以及genes.gff，config文件添加物種命名

mkdir snpEFF

mkdir data

mkdir Cucumber (對黃瓜進行注釋，專門地文件夾)

拷貝基因組文件（復(fù)制為sequences.fa）

cp /share/home/bmk366/reseq/test-data/genome/cucumber_ChineseLong_v2_genome.fa sequences.fa

拷貝gff文件（復(fù)制為genes.gff）

cp /share/home/bmk366/reseq/test-data/genome/cucumber_ChineseLong_v2.gff3 genes.gff


vim snpEff.config

在文件結(jié)尾添加一行內(nèi)容

Cucumber.genome: Cucumber（添加內(nèi)容與文件名相一致）

##構(gòu)建物種數(shù)據(jù)庫

java -jar /share/nas2/genome/biosoft/snpEff/snpEff.jar build -c ./snpEff.config -gff3 Cucumber

###進行注釋

java -jar /share/nas2/genome/biosoft/snpEff/snpEff.jar eff Cucumber ../variant_calling/raw.filter.snp.vcf -c ./snpEff.config -o gatk -ud 5000 >./raw.snp.anno.vcf

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

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

  • ### 一般從測序公司拿到的下級數(shù)據(jù)有raw data，或者經(jīng)過質(zhì)控之后的clean data，格式為fq。參考基...
    纖蹄馬閱讀 898評論 0贊 7
  • 第三節(jié)主要內(nèi)容：測序?qū)嶒灹鞒?、測序原理及基本名詞解釋 1. 測序錯誤率原因：Phasing & Pre-phasi...
    不想透明的小透明閱讀 3,902評論 0贊 5
  • 通過雙端測序的數(shù)據(jù)比對參考基因組。 第一步準(zhǔn)備工作（這里只用了擬南芥1號染色體的數(shù)據(jù)） bwa index ./s...
    每天都想睡覺的阿源閱讀 1,938評論 0贊 6
  • 我通過查資料獲得已知達松維爾擬諾卡氏菌亞種（cardiopsis dassonvillei subsp. dass...
    lizg閱讀 11,848評論 2贊 37
  • 基因組重測序數(shù)據(jù)目的：需要檢測基因組中的變異，找到并定位這些突變位點 條件：參考基因組、重測序數(shù)據(jù)、 分析流程： ...
    古月福_閱讀 3,121評論 1贊 2

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