2021-01-06GWAS全基因组关联分析流程


转载自:https://blog.csdn.net/yangl7/article/details/109202042

一、BWA比对

1.构建索引

bwa index -a is example.fasta #构建索引 -a is算法 (BWT构造算法:bwtsw、is或rb2)

2.进行比对

bwa mem -t 6 -R '@RG\tID:foo\tPL:Illumina\tSM:example'

example.fasta example_1.fq.gz example_2.fq.gz > example.sam

# 进行对比 mem算法 -t 运行的核数目

# -R添加头部 

ID:这是Read Group的分组ID,一般设置为测序的lane ID(不同lane之间的测序过程认为是独立的),下机数据中我们都能看到这个信息的,一般都是包含在fastq的文件名中;

PL:指的是所用的测序平台,这个信息不要随便写,在GATK中,PL只允许被设置为:ILLUMINA,SLX,SOLEXA,SOLID,454,LS454,COMPLETE,PACBIO,IONTORRENT,CAPILLARY,HELICOS或UNKNOWN这几个信息。如果不知道,那么必须设置为UNKNOWN。

SM:样本ID。

LB:测序文库的名字,如果上面的lane ID足够用于区分的话,也可以不用设置LB;

(用GATK检测变异 其中ID,PL和SM信息是必须的)

二、samtools格式转换

1.sam格式转换为bam格式

samtools view -bS example.sam -o example.bam

# -b 输出bam格式文件 -S 输入sam格式文件

2.质控

samtools view -h -b -q30 example.bam > example.q30.bam

# -q 比对的最低质量值 -h 输出的文件包含头部信息 -b 输出bam格式文件

3.构建索引

samtools faidx base/example.fasta

# 该命令会在example.fasta所在目录下创建一个example.fai索引文件

gatk CreateSequenceDictionary -R example.fasta -O example.dict

# 创建gatk索引 生产dict文件

三、gatk变异检测

1.排序

gatk SortSam -I example.q30.bam -O example.q30.sort.bam

-R base/example.fasta -SO coordinate --CREATE_INDEX true

# -I 输入文件 -O 输出文件 -R参考基因组 --CREATE_INDEX 是否建立索引

将sam文件中同一染色体对应的条目按照坐标顺序从小到大进行排序

2.标记重复序列

gatk MarkDuplicates -I example.q30.sort.bam -O example.q30.sort.markdup.bam -M example.q30.sort.markdup_metrics.txt

# -I 输入排序后的文件 -O 输出文件 -M 输出重复矩阵

注意:

samtools index example.q30.sort.markdup.bam

# 为gatk建立索引这里非常重要

3.检测变异

gatk HaplotypeCaller --emitRefConfidence GVCF -R base/example.fasta -I example.q30.sort.markdup.bam -L chrX -O chrX.g.vcf.gz

# HaplotypeCaller同时检测snp和indel -R 参考基因组 -I 输入文件 -L 仅检测该染色体的变异(分染色体检测变异,加快速度)-O 输出文件

这里分染色体进行检测,后续再进行合并,可以加快检测速度。

合并文件(g.vcf)

gatk CombineGVCFs -R base/example.fasta --variant example1.g.vcf.gz --variant example2.g.vcf.gz -O con.vcf.gz

# -R 参考基因组 --variant 输入变异文件 可以输入多个文件 -O 输出文件

检测变异

gatk GenotypeGVCFs -R ref.fa -V test.g.vcf -O test.vcf

4.提取SNP变异

gatk SelectVariants -R base/example.fasta -V test.vcf -O test.snp.vcf --select-type-to-include SNP

# -R 参考基因组 -O 输出vcf文件 -V 输入vcf文件 --select-type-to-include 选取提取的变异类型(#SNP,MNP,INDEL,SYMBOLIC,MIXED)

5.对SNP进行过滤

gatk VariantFiltration -O chr5.Filt.vcf -V chr5.vcf

--cluster-window-size 10

--filter-expression "MQ0 >= 4 && ((MQ0 / (1.0 * DP)) > 0.1)" --filter-name "HARD_TO_VALIDATE"

--filter-expression "DP < 5 " --filter-name "LowCoverage"

--filter-expression "QUAL < 30.0 " --filter-name "VeryLowQual"

--filter-expression "QUAL > 30.0 && QUAL < 50.0 "

--filter-name "LowQual" --filter-expression "QD < 1.5 "

--filter-name "LowQD"

# -O 输出文件 -V输入变异文件 --filter-expression 过滤条件 --filter-name 被过滤掉的SNP不会删除,而是给一个标签,例如 "LowCoverage" 。

这里可以将过滤条件合并,仅给出一个标签。--cluster-window-size 以10个碱基为一个窗口

这里通过设定相应的参数值进行了硬过滤,实际应用时还要根据数据类型及自己的需求设定相应的参数。

6.合并文件(vcf)

删除掉被过滤的SNP

grep -v "LowCoverage" Filt.vcf > Filt1.vcf

# -v显示不包含匹配文本的所有行 "LowCoverage"上一步给出的标签

合并成一个文件

gatk MergeVcfs -I chr1.Filt.vcf -I chr2.Filt.vcf -I chr3.Filt.vcf  -O Filt.vcf

# -I 输入文件 可以多次指定 -O 输出文件

至此为止就得到了整个群体的VCF变异文件,后续都是基于此文件来进行相应的分析。

四、Plink格式转换及主成分分析

1.VCF格式转换为 ped/map格式

vcftools --vcf snp.vcf --plink --out snp

2.ped/map格式转换为bed/bim/fam格式

plink --file snp --make-bed --out snp_test

3.主成分分析

Plink --threads 8 --bfile snp --pca 10 --out pca

# --threads 线程数 --bfile 输入.bed文件 --pca 主成分的成分数 --out输出的文件名

五、Admixture 群体结构

1.群体结构分析

for K in 2 3 4 5 6 7 8 9 10; \

do admixture --cv hapmap3.bed $K | tee log${K}.out; done

#2 3 4 5 6 7 8 9 10分成的群体结构数 hapmap3.bed 输入文件

注意:

如果你的数据格式是plink的bed文件, 比如a.bed, 那么你应该包含a.bim, a.fam

如果你的数据格式是plink的ped文件, 比如b.ped, 那么你应该包括b.map

K值根据实际情况进行设置,通过比较得到最佳K值。

grep -h CV log*.out

#查看最佳K值 输出最佳K值文件:hapmap3.3.Q

2.R语言作图

tbl=read.table("hapmap3.3.Q")

barplot(t(as.matrix(tbl)), col=rainbow(K),

xlab="Individual", ylab="Ancestry", border=NA)

hapmap.3.Q文件为群体结构的结果,作为协变量进行关联分析

六、Tassel关联分析

Tassel的管道命令不允许有回车符号,使用以下命令时需要将#注释及换行删除。

1.VCF格式转换为 hmp格式

run_pipeline.pl -SortGenotypeFilePlugin -inputFile example.vcf -outputFile example -fileType VCF

#给vcf文件排序,排成tassel认可的序列

#-inputFile 输入的文件名 -outputFile 输出的文件名 -fileType 输出的文件格式

run_pipeline.pl -fork1 -vcf example.vcf  -export example -exportType Hapmap -runfork1

#vcf文件转换为hapmap格式

#-vcf 输出的文件 -export 输出的文件 -exportType 输出的文件格式

2.亲缘关系

run_pipeline.pl -importGuess genotype.hmp.txt #打开数据文件

-KinshipPlugin -method Centered_IBS -endPlugin #计算亲缘关系

-export genotype_kinship.txt  #输出文件名

-exportType SqrMatrix #输出格式

3.关联分析

混合线性模型

run_pipeline.pl -fork1 -h genotype.hmp.txt -filterAlign -filterAlignMinCount 150 -filterAlignMinFreq 0.05 -filterAlignMaxFreq 1

#-h读取hapmap格式的基因型数据 -filterAlign 打开过滤选项

#-filterAlignMinCount -filterAlignMinFreq -filterAlignMaxFreq 过滤的条件 需要按照实际情况进行修改

-fork2 -r traits.txt

#-r 读取表型数据

-fork3 -q population_structure.txt -excludeLastTrait

#-q 读取群体结构数据 -excludeLastTrait 去掉最后一个群体结构 (当分成的群体结构>2时)

-fork4 -k kinship.txt

#读取亲缘关系数据

-combine5 -input1 -input2 -input3 -intersect -combine6 -input5 -input4 -mlm -export example

#-combine合并数据 -mlm混合线性模型 -export输出文件名

一般线性模型

run_pipeline.pl

-fork1 -h genotype.hmp.txt -filterAlign

-filterAlignMinCount 150 -filterAlignMinFreq 0.05 -filterAlignMaxFreq 1

-fork2 -r traits.txt

-fork3 -q population_structure.txt -excludeLastTrait

-combine4 -input1 -input2 -input3 -intersect -glm -export example

#-glm 一般线性模型

4.R语言作图

Library(qqman)

#加载qqman包

曼哈顿图

manhattan(example,ylim=c(0,10),col = color_set,annotatePval = 0.01)

# ylim Y轴范围 col 颜色 annotatePval 标记最高位点 CHR==1 绘制每个染色体的曼哈顿图

Q-Q plot

qq(example$P)

七、其他

1.基因组统计工具

可以统计fasta和fastq文件中的信息。

seqkitfx2tab example.fasta -l -n

-l统计序列长度 -n 统计染色体

2.提取文本文档中某列

用于Tassel关联分析后的结果文件,提取相应的列进行R语言绘图。

catMLM.txt | awk '{print $1" "$3" "$4" "$7}' > manhattan.txt

# $提取的列数

3.删除文本文档中不包含匹配文本的行

用于过滤后删除低质量的SNP。

grep-v"LowCoverage"Filt.vcf > Filt1.vcf

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 206,378评论 6 481
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 88,356评论 2 382
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 152,702评论 0 342
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 55,259评论 1 279
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 64,263评论 5 371
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,036评论 1 285
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,349评论 3 400
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,979评论 0 259
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 43,469评论 1 300
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,938评论 2 323
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,059评论 1 333
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,703评论 4 323
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,257评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,262评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,485评论 1 262
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,501评论 2 354
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,792评论 2 345

推荐阅读更多精彩内容