美文网首页GWASQTL定位R语言可视化之美
Post-GWAS-LD图与基因,关联结果一起可视化

Post-GWAS-LD图与基因,关联结果一起可视化

作者: Hello育种 | 来源:发表于2021-10-28 05:21 被阅读0次
    image.png

    IntAssoplot包可以将GWAS结果和LD, 基因结构图整合在一起画出。
    介绍来自:https://github.com/whweve/IntAssoPlot

    1下载IntAssoplot包

    #### 前期准备包
    install.packages(c("devtools","remotes"))
    
    #install depended packages, including ggplot2, SNPRelate, ggrepel, gdsfmt and reshape2 #ggplot2, ggrepel, and reshape2 are installed from CRAN
    
    install.packages(c("ggplot2","ggrepel","reshape2"))
    
    #SNPRelate and gdsfmt are installed from Bioconductor
    
    if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
    
    BiocManager::install(c("SNPRelate","gdsfmt"))
    ##### 下载IntAssoPlot from Github
    library(remotes) # version 2.1.0
    
    #download, build, and install IntAssoPlot without creating vignette
    
    install_github("whweve/IntAssoPlot")
    或者
    #download, build, and install IntAssoPlot with creating vignette
    
    install_github("whweve/IntAssoPlot",build=TRUE,build_vignettes = TRUE)
    

    2 数据准备

    2.1 GWAS的结果文件

    包含关联分析文件的数据集。 目前,我们针对相应标记的物理位置绘制了 -log10 转换后的 p 值,该值可以从基因组扫描中得出。 因此,需要一个包含标记名称及其 p 值的数据框。 需要的变量如下:Marker(分子标记名称)、Locus(标记的染色体)、Site(标记的位置)、p(标记的p值)这里是关联数据格式的快速演示:

    head(association)
    #>     Marker  Locus     Site         p
    #> 1 SNP37049     9 90042244 0.2073684
    #> 2 SNP37050     9 90042408 0.5028760
    #> 3 SNP37052     9 90210602 0.1561271
    #> 4 SNP37053     9 90255652 0.2814653
    #> 5 SNP37056     9 90290939 0.6392511
    #> 6 SNP44597     9 90367887 0.6582597
    #the attribute of each column could be viewed as:
    str(association$Marker)
    #>  chr [1:2316] "SNP37049" "SNP37050" "SNP37052" "SNP37053" "SNP37056" ...
    str(association$Locus)
    #>  int [1:2316] 9 9 9 9 9 9 9 9 9 9 ...
    str(association$Site)
    #>  int [1:2316] 90042244 90042408 90210602 90255652 90290939 90367887 90377686 90385837 90560965 90759605 ...
    str(association$p)
    #>  num [1:2316] 0.207 0.503 0.156 0.281 0.639 ..
    

    2.2 gene结构注释的文件格式

    包含注释文件的数据集,通常是 gtf 文件,没有列名。 当对一个物种的基因组进行测序时,可以在基因组注释网站上找到 gtf 文件。 对于已发布的基因组,请参阅 www.ensembl.org。 由于大多数 gtf 文件缺少列名,因此可以使用 read.table(?where is your file?,header=FALSE) 读取注释文件。 这是 gtf 数据格式的快速演示:

    head(gtf)
    #>   V1             V2          V3       V4       V5 V6 V7 V8
    #> 1  9 protein_coding        exon 90030791 90030847  .  -  .
    #> 2  9 protein_coding        exon 90030635 90030701  .  -  .
    #> 3  9 protein_coding         CDS 90030635 90030681  .  -  0
    #> 4  9 protein_coding start_codon 90030679 90030681  .  -  0
    #> 5  9 protein_coding        exon 90029368 90029471  .  -  .
    #> 6  9 protein_coding         CDS 90029368 90029471  .  -  1
    #>                                                                                                             V9
    #> 1                gene_id "GRMZM2G416644"; transcript_id "GRMZM2G416644_T01"; exon_number "1"; seqedit "false";
    #> 2                gene_id "GRMZM2G416644"; transcript_id "GRMZM2G416644_T01"; exon_number "2"; seqedit "false";
    #> 3 gene_id "GRMZM2G416644"; transcript_id "GRMZM2G416644_T01"; exon_number "2"; protein_id "GRMZM2G416644_P01";
    #> 4                                 gene_id "GRMZM2G416644"; transcript_id "GRMZM2G416644_T01"; exon_number "2";
    #> 5                gene_id "GRMZM2G416644"; transcript_id "GRMZM2G416644_T01"; exon_number "3"; seqedit "false";
    #> 6 gene_id "GRMZM2G416644"; transcript_id "GRMZM2G416644_T01"; exon_number "3"; protein_id "GRMZM2G416644_P01";
    #the attribute of each column could be viewed as:
    str(gtf$V1)
    #>  chr [1:4665] "9" "9" "9" "9" "9" "9" "9" "9" "9" "9" "9" "9" "9" "9" "9" ...
    str(gtf$V2)
    #>  chr [1:4665] "protein_coding" "protein_coding" "protein_coding" ...
    str(gtf$V3)
    #>  chr [1:4665] "exon" "exon" "CDS" "start_codon" "exon" "CDS" "exon" "CDS" ...
    str(gtf$V4)
    #>  int [1:4665] 90030791 90030635 90030635 90030679 90029368 90029368 90028806 90029227 90029224 90030697 ...
    str(gtf$V5)
    #>  int [1:4665] 90030847 90030701 90030681 90030681 90029471 90029471 90029282 90029282 90029226 90031013 ...
    str(gtf$V6)
    #>  chr [1:4665] "." "." "." "." "." "." "." "." "." "." "." "." "." "." "." ...
    str(gtf$V7)
    #>  chr [1:4665] "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" ...
    str(gtf$V8)
    #>  chr [1:4665] "." "." "0" "0" "." "1" "." "2" "0" "." "0" "0" "." "0" "0" ...
    str(gtf$V9)
    #>  chr [1:4665] "gene_id \"GRMZM2G416644\"; transcript_id \"GRMZM2G416644_T01\"; exon_number \"1\"; seqedit \"false\";" ...
    

    2.3 基因型数据的格式

    包含基因型文件(通常是 hapmap 文件)和列名的数据集。 Hapmap 是一种常用的存储基因信息的格式。 可以对重要材料进行高通量测序,将序列与参考基因组对齐,提取 SNP/InDels。 另一方面,人们可以对一个特定的基因进行重新测序。 为此,设计重叠扩增引物、扩增基因组或转录组片段、多重比对序列、提取 SNP/InDels。 这是关联数据格式的快速演示:

    #only 20 column of the genotype markers are shown.
    head(zmvpp1_hapmap[,1:20])
    #>          rs       allele   chrom  pos strand assembly center protLSID assayLSID
    #> 1 INDEL-665    -/+     9 94183611      +       NA     NA       NA        NA
    #> 2   snp-663    A/T     9 94183597      +       NA     NA       NA        NA
    #> 3   snp-649    A/G     9 94183594      +       NA     NA       NA        NA
    #> 4   snp-646    C/T     9 94183585      +       NA     NA       NA        NA
    #> 5   snp-637    G/T     9 94183584      +       NA     NA       NA        NA
    #> 6   snp-636    A/G     9 94183575      +       NA     NA       NA        NA
    #>   panel QCcode B73 CIMBL55 CIMBL32 CML298 CML170 CIMBL12 CIMBL95 CIMBL1 CIMBL56
    #> 1    NA     NA  --      ++      ++     ++     ++      ++      ++     --      --
    #> 2    NA     NA  AA      TT      TT     TT     TT      TT      TT     NN      NN
    #> 3    NA     NA  AA      GG      GG     GG     GG      GG      GG     NN      NN
    #> 4    NA     NA  CC      TT      TT     TT     TT      TT      TT     NN      NN
    #> 5    NA     NA  GG      TT      TT     TT     TT      TT      TT     NN      NN
    #> 6    NA     NA  AA      AA      AA     AA     AA      AA      AA     NN      NN
    

    3 绘制与注释和 LD 矩阵的关联

    IntAssoPlot 是将关联分析、基因结构和 LD 矩阵的结果自动整合到一个视图中。 事实上,这项任务很困难,因为散点图 - log10P 与 SNP 的物理位置相对,而 LD 矩阵没有物理位置。该包的工作其实就是为了解决上述问题。 在这里,我们使用先前发布的数据(Wang 等人,2016 年)展示了一个示例来展示 IntAssoPlot 的用法。

    3.1 具有一组基因型标记的区域综合图

    在跨越 400 kbp 区域的区域绘制关联结果,并使用与关联映射相同的 SNP 标记绘制 LD 矩阵。

    IntRegionalPlot(chr=9,left=94178074-200000,right=94178074+200000,gtf=gtf,association=association,hapmap=hapmap_am368,hapmap_ld=hapmap_am368,threshold=5,leadsnp_size=2)
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 368
    #>     # of SNPs: 270
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 143276
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 368
    #>     # of SNPs: 270
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 143276
    

    结果图:


    image.png

    3.2 用从浅灰色到深灰色的颜色绘制 LD 值

    IntRegionalPlot(chr=9,left=94178074-200000,right=94178074+200000,gtf=gtf,association=association,hapmap=hapmap_am368,hapmap_ld=hapmap_am368,threshold=5,leadsnp_size=2,colour02 = "gray1",colour04 = "gray21",colour06 = "gray41",colour08 = "gray61",colour10 = "gray81",)
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 368
    #>     # of SNPs: 270
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 143276
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 368
    #>     # of SNPs: 270
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 143276
    
    image.png

    3.3 用从白色到红色的颜色绘制 LD 值

    #get five colors ranging from white to red
    pal <- colorRampPalette(c("white", "red"))
    IntRegionalPlot(chr=9,left=94178074-200000,right=94178074+200000,gtf=gtf,association=association,hapmap=hapmap_am368,hapmap_ld=hapmap_am368,threshold=5,leadsnp_size=2,colour02 = pal(5)[1],colour04 = pal(5)[2],colour06 = pal(5)[3],colour08 = pal(5)[4],colour10 = pal(5)[5])
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 368
    #>     # of SNPs: 270
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 143276
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 368
    #>     # of SNPs: 270
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 143276
    
    image.png

    3.4 用从白色到红色的颜色绘制 LD 值并标记基因名称

    #get five colors ranging from white to red
    pal <- colorRampPalette(c("white", "red"))
    IntRegionalPlot(chr=9,left=94178074-200000,right=94178074+200000,gtf=gtf,association=association,hapmap=hapmap_am368,hapmap_ld=hapmap_am368,threshold=5,leadsnp_size=2,colour02 = pal(5)[1],colour04 = pal(5)[2],colour06 = pal(5)[3],colour08 = pal(5)[4],colour10 = pal(5)[5],label_gene_name = TRUE)
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 368
    #>     # of SNPs: 270
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 143276
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 368
    #>     # of SNPs: 270
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 143276
    
    image.png

    3.5 具有两组基因型标记的区域综合图

    在跨越 200 kbp 区域的区域绘制关联结果,并使用与关联映射不同的 SNP 标记器绘制 LD 矩阵。 此功能允许研究人员在更广泛的标记范围内研究 LD 结构。

    IntRegionalPlot(chr=9,left=94178074-100000,right=94178074+100000,gtf=gtf,association=association,hapmap=hapmap_am368,hapmap_ld=hapmap2,threshold=5,leadsnp_size=2)
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 368
    #>     # of SNPs: 108
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 56746
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 104
    #>     # of SNPs: 1,081
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 114672
    
    image.png

    3.6 具有一组基因型标记的相对较小的区域综合图

    IntRegionalPlot(chr=9,left=94178074-2000,right=94178074+5000,gtf=gtf,association=association,hapmap=hapmap_am368,hapmap_ld=hapmap_am368,threshold=5,leadsnp_size=2)
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 368
    #>     # of SNPs: 41
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 21412
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 368
    #>     # of SNPs: 41
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 21412
    
    image.png

    3.7 单基因水平图

    绘制给定基因的关联结果,并使用与关联映射相同的 SNP 标记物绘制 LD 矩阵。 还指定的标记通过各种形状和颜色突出显示。

    3.7.1 基本版

    IntGenicPlot('GRMZM2G170927_T01',gtf,association=zmvpp1_association,hapmap=zmvpp1_hapmap,hapmap_ld = zmvpp1_hapmap,threshold=8,leadsnpLD = FALSE)
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 141
    #>     # of SNPs: 53
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 9992
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 141
    #>     # of SNPs: 53
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 9992
    
    image.png

    3.7.2从基因的上游/下游扩展区域

    IntGenicPlot('GRMZM2G170927_T01',gtf,association=zmvpp1_association,hapmap=zmvpp1_hapmap,hapmap_ld = zmvpp1_hapmap,threshold=8,up=500,down=600,leadsnpLD = FALSE)
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 141
    #>     # of SNPs: 124
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 23488
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 141
    #>     # of SNPs: 124
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 23488
    
    image.png

    3.7.3 突出显示选定的标记,在数据框中指定颜色和形状:标记 2 突出显示

    IntGenicPlot('GRMZM2G170927_T01',gtf,association=zmvpp1_association,hapmap=zmvpp1_hapmap,hapmap_ld = zmvpp1_hapmap,threshold=8,up=500,down=600,leadsnpLD = FALSE,marker2highlight=marker2highlight)
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 141
    #>     # of SNPs: 124
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 23488
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 141
    #>     # of SNPs: 124
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 23488
    
    image.png

    3.7.4 添加连接线

    IntGenicPlot('GRMZM2G170927_T01',gtf,association=zmvpp1_association,hapmap=zmvpp1_hapmap,hapmap_ld = zmvpp1_hapmap,threshold=8,up=500,down=600,leadsnpLD = FALSE,marker2highlight=marker2highlight,link2gene=marker2link,link2LD=marker2link)
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 141
    #>     # of SNPs: 124
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 23488
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 141
    #>     # of SNPs: 124
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 23488
    
    image.png

    3.7.5 为高亮标记添加名称

    IntGenicPlot('GRMZM2G170927_T01',gtf,association=zmvpp1_association,hapmap=zmvpp1_hapmap,hapmap_ld = zmvpp1_hapmap,threshold=8,up=500,down=600,leadsnpLD = FALSE,marker2highlight=marker2highlight,link2gene=marker2link,link2LD=marker2link,marker2label=marker2link,marker2label_angle=60,marker2label_size=2)
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 141
    #>     # of SNPs: 124
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 23488
    #> Linkage Disequilibrium (LD) estimation on genotypes:
    #>     # of samples: 141
    #>     # of SNPs: 124
    #>     using 1 thread
    #>     method: correlation
    #> LD matrix:    the sum of all selected genotypes (0,1,2) = 23488
    
    image.png

    相关文章

      网友评论

        本文标题:Post-GWAS-LD图与基因,关联结果一起可视化

        本文链接:https://www.haomeiwen.com/subject/pqmmaltx.html