GWAS-3K分析

  1. 将现有名称进行统一化处理

    第一种办法是从网站下载所有品种的名称信息，从中筛选目标470份样品的名称I找到对应的IRIS名称

    第二种办法是利用my_list菜单提交自己的IRIS号，或其他的信息

  2. 利用plink提取目标品种的SNP信息

    当前目录下，存在bed/bim/fam三个文件

    plink --bfile pruned_v2.1 --noweb --keep iris_id --recode --make-bed --out cold

iris_id含有两列，一列为fid，一列为iid，此处相同

plink --bfile cold --export vcf-iid --out cold

将bed/bim/fam生成vcf文件

plink --vcf in.vcf --maf 0.05 --geno 0.2 --mind 0.1 --recode vcf-iid --out in.filter --allow-extra-chr

-iid生成的文件不含fid，也不含有下划线

conda 安装tassel,生成hmp格式的文件

nohup ... & 后台执行

nohup run_pipeline.pl -Xms10g -Xmx100g -vcf cold.filter.vcf -sortPositions -export hmp.txt -exportType HapmapDiploid &

利用admixture进行群体结构分析,最小的CV值对应的k值较好，时间较长

#PBS -N admixture
#PBS -l nodes=1:ppn=8
#PBS -q middle

cd $PBS_O_WORKDIR

for K in {1..10}; do admixture --cv  cold.filter.bed $K | tee log${K}.out; done

grep -h CV log*.out
#群体结构绘图，选择最小cv对应的Q矩阵，以及可视化网站
#http://omicsspeaks.com/strplot2/

paste <(cut -f1 demo.admixture.in.nosex) demo.admixture.in.9.Q |sed 's/^/G\t/g' |sed 		's/\s/,/g'>plotQ.csv

#服务器上运行R脚本
R --slave --vanilla <mrmlm.r
R --slave --vanilla <rMVP.R
#--slave 只打印输出，不打印过程--vanilla整合多种参数

​

rMVP包进行关联分析，对应模型为GLM，MLM，FarmCPU

github说明

library('rMVP')
MVP.Data(filePhe='phenotype.txt',fileVCF='cold.filter.vcf',out='mytest')
dir()
#generates 5 files
#load the data 
pheno<-read.table('mytest.phe',header=TRUE)
geno<-attach.big.matrix('mytest.geno.desc')
map<-read.table('mytest.geno.map',header=TRUE)
#if you have more than one phentype,use follow method
for(i in 2:ncol(pheno)){
  imMVP <- MVP(
    phe=pheno[, c(1, i)],
    geno=geno,
    map=map,
    #K=Kinship,
    #CV.GLM=Covariates,
    #CV.MLM=Covariates,
    #CV.FarmCPU=Covariates,
    nPC.GLM=5,
    nPC.MLM=3,
    nPC.FarmCPU=3,
    priority="speed",
    #ncpus=10,
    vc.method="BRENT",
    maxLoop=10,
    method.bin="static",
    #permutation.threshold=TRUE,
    #permutation.rep=100,
    threshold=0.05,
    method=c("GLM", "MLM", "FarmCPU")
  )
  gc()
}

mrMLM分析

#linux environment
# remember chmod 755 mrmlm.R

#!/usr/local/bin/Rscript
library('mrMLM')
mrMLM(fileGen='hmp.txt',filePhe='phenotype.txt',Genformat='Hmp',method=c("mrMLM","FASTmrMLM","FASTmrEMMA","pLARmEB","pKWmEB","ISIS EM-BLASSO"),trait=1:12,CriLOD=3,DrawPlot=TRUE,Plotformat='jpeg',dir='/home/caisl/lee/gwas/cold/mrMLM')

LD分析

~/biosoft/LDBlockShow/bin/LDBlockShow -InVCF sample.filter.vcf -OutPut out -Region 2:2000-10000 -OutPng -SeleVar 1
#若染色体名称为chr2,就改为chr2：2000-10000

群体内个体间遗传距离分析

plink --file hapmap --cluster --matrix --noweb
#默认情况下会生成plink.mibs文件，是一个距离矩阵，可以用R语言读取，然后进行MDS分析。

plink --file hapmap --cluster --mds-plot 10
#Multidimensional scaling plots（MDS） 实际上是多维变量分解IBS矩阵。这里的MDS和PCA是一样的。


#绘图部分,可以直接利用ibs距离进行分析
m <- as.matrix(read.table("plink.mibs"))
mds <- cmdscale(as.dist(1-m))
k <- c( rep("green",45) , rep("blue",44) )
plot(mds,pch=20,col=k)
legend("topleft", legend = c("CHB", "JPB"), pch = 20, col = c("green", "blue"))
#利用mds的前两个主成分（维度）分析
library('ggplot2')
library(ggrepel)
dat <- read.table ( "plink.mds", header = T,stringAsFactors=FALSE)
data<-dat[,c(1,4,5)]
C0<-data$FID
C1<-data$C1
C2<-data$C2
head(data)
ggplot(data)+geom_point(aes(C1,C2),color='red')+geom_text_repel(aes(data$C1,data$C2,label=data$FID))+labs(main = 'MDS PLOT')


#运行GAPIT包

options('repos')#查看当前镜像
options(repos="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")#修改镜像
source("http://zzlab.net/GAPIT/GAPIT.library.R")
#载入需要的程辑包：ape
#载入需要的程辑包：compiler
#载入需要的程辑包：EMMREML
#载入需要的程辑包：Matrix
#3载入需要的程辑包：scatterplot3d
myY<-read.table('468phenotype.txt',header=TRUE, stringsAsFactors=FALSE)
myG<-read.delim('hmp.txt',header=FALSE, stringsAsFactors=FALSE)
source("http://zzlab.net/GAPIT/gapit_functions.txt")
mtGAPIT<-GAPIT(Y=myY,G=myG,PCA.total=2,model=c('CMLM','Blink'))
#在自动生成主成分图片和kinship矩阵