重測序分析

重測序：是對已對已知基因組的物種進(jìn)行測序枫疆，去挖掘不同個體和群體之間的差異性闷沥。

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

SNP，INDEL, SV , SNV

進(jìn)化分析善茎，群體結(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

實際進(jìn)行比對

/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文件按照比對到參考基因組的順序進(jìn)行排序：

/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擴增或者錯誤測序贞绳，進(jìn)行標(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ā)生錯配的地方進(jìn)行重新比對提高假陽性和假陰性萎攒。

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進(jìn)行劃窗分析

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 位點進(jìn)行過濾

-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位點的文件

###使用SNPEFF 對SNP 位點進(jìn)行注釋

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

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

mkdir snpEFF

mkdir data

mkdir Cucumber (對黃瓜進(jìn)行注釋束亏，專門地文件夾)

拷貝基因組文件（復(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

###進(jìn)行注釋

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)系作者

人面猴
序言：七十年代末，一起剝皮案震驚了整個濱河市供汛，隨后出現(xiàn)的幾起案子枪汪，更是在濱河造成了極大的恐慌涌穆，老刑警劉巖怔昨，帶你破解...
沈念sama閱讀 206,013評論 6贊 481
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件，死亡現(xiàn)場離奇詭異宿稀，居然都是意外死亡趁舀，警方通過查閱死者的電腦和手機，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 88,205評論 2贊 382
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進(jìn)店門祝沸，熙熙樓的掌柜王于貴愁眉苦臉地迎上來矮烹，“玉大人，你說我怎么就攤上這事罩锐》畋罚” “怎么了？”我有些...
開封第一講書人閱讀 152,370評論 0贊 342
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵涩惑，是天一觀的道長仁期。經(jīng)常有香客問我，道長竭恬，這世上最難降的妖魔是什么跛蛋？我笑而不...
開封第一講書人閱讀 55,168評論 1贊 278
?港島之戀（遺憾婚禮）
正文為了忘掉前任，我火速辦了婚禮痊硕，結(jié)果婚禮上赊级，老公的妹妹穿的比我還像新娘。我一直安慰自己岔绸，他們只是感情好理逊，可當(dāng)我...
茶點故事閱讀 64,153評論 5贊 371
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布橡伞。她就那樣靜靜地躺著，像睡著了一般晋被。火紅的嫁衣襯著肌膚如雪骑歹。梳的紋絲不亂的頭發(fā)上，一...
開封第一講書人閱讀 48,954評論 1贊 283
城市分裂傳說
那天墨微，我揣著相機與錄音道媚，去河邊找鬼。笑死翘县，一個胖子當(dāng)著我的面吹牛最域，可吹牛的內(nèi)容都是我干的。我是一名探鬼主播锈麸，決...
沈念sama閱讀 38,271評論 3贊 399
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼镀脂，長吁一口氣：“原來是場噩夢啊……” “哼！你這毒婦竟也來了忘伞？” 一聲冷哼從身側(cè)響起薄翅，我...
開封第一講書人閱讀 36,916評論 0贊 259
萬榮殺人案實錄
序言：老撾萬榮一對情侶失蹤，失蹤者是張志新（化名）和其女友劉穎氓奈，沒想到半個月后翘魄，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體，經(jīng)...
沈念sama閱讀 43,382評論 1贊 300
?護(hù)林員之死
正文獨居荒郊野嶺守林人離奇死亡舀奶，尸身上長有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點故事閱讀 35,877評論 2贊 323
?白月光啟示錄
正文我和宋清朗相戀三年暑竟，在試婚紗的時候發(fā)現(xiàn)自己被綠了。大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片育勺。...
茶點故事閱讀 37,989評論 1贊 333
活死人
序言：一個原本活蹦亂跳的男人離奇死亡但荤，死狀恐怖，靈堂內(nèi)的尸體忽然破棺而出涧至，到底是詐尸還是另有隱情腹躁，我是刑警寧澤，帶...
沈念sama閱讀 33,624評論 4贊 322
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布南蓬，位于F島的核電站纺非，受9級特大地震影響，放射性物質(zhì)發(fā)生泄漏蓖康。R本人自食惡果不足惜铐炫，卻給世界環(huán)境...
茶點故事閱讀 39,209評論 3贊 307
男人毒藥：我在死后第九天來索命
文/蒙蒙一、第九天我趴在偏房一處隱蔽的房頂上張望蒜焊。院中可真熱鬧倒信，春花似錦、人聲如沸泳梆。這莊子的主人今日做“春日...
開封第一講書人閱讀 30,199評論 0贊 19
一樁弒父案榜掌，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽。三九已至乘综，卻和暖如春憎账，著一層夾襖步出監(jiān)牢的瞬間，已是汗流浹背卡辰。一陣腳步聲響...
開封第一講書人閱讀 31,418評論 1贊 260
情欲美人皮
我被黑心中介騙來泰國打工胞皱，沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留，地道東北人九妈。一個月前我還...
沈念sama閱讀 45,401評論 2贊 352
代替公主和親
正文我出身青樓反砌，卻偏偏與公主長得像，于是被迫代替她去往敵國和親萌朱。傳聞我的和親對象是個殘疾皇子宴树，可洞房花燭夜當(dāng)晚...
茶點故事閱讀 42,700評論 2贊 345

重測序分析

推薦閱讀更多精彩內(nèi)容