SNP分析命令

E:\ > cd e:

E:\

E:\ > cd plink-1

E:\plink-1>plink –file test

1.Map 更新

Plink --sheep --file data --update-map position.txt --recode --out data1

Chrnew.txt -- update-chr --recode --out data2 Position: SNP code and position

Chrnew：SNP code and Chr.

2．SNP merge

Plink --file data1 --merge data2.ped data2.map --recode --out merge

3．提取SNP位点

Plink --file data --extract 50kSNP.txt --recode --out data1

50kSNP.txt: 50k中的SNP名

4. Quality control

Call rate >98%/99%

Plink --file sheep --geno 0.02 --recode --out sheepgeno

Plink --file sheepgeno --mind 0.01 --recode --out sheepmind

MAF>0.05

Plink --file sheepmind --maf 0.05 --recode --out sheepmaf

Hardy-Weinberg equilibrium <0.0001

Plink --file sheepmaf --hwe 0.0001 --recode --out sheephwe

Exclude the SNP markers with either chromosome or both unknown

Plink --sheep --file sheephwe --extract 4newsnp.txt --recode --out sheep4

Note: 制作4newsnp.txt（包含chromosome 和base-pair position 都为0的SNP）

To identify sample duplication or half-sibs or closer

Plink –sheep –file sheep4 –genome –max 0.85

Note：Check the genome file

5. LD quality control

Plink –sheep --file sheep4 –indep-pairwise 100 25 0.2 –out sheepld0.2

Plink --sheep --file sheep4 --indep-pairwise 100 25 0.05 --out sheepld0.05

Plink--file sheep4--ld-window-r2 0.2 --out sheepldr0.2

输出结果为data prunein 和data prune out

(质控时，要去除X染色体)

将data prune in 转化为ped和map

Plink --sheep --file 114hwe --extract 114sheep0.05.prune.in --recode --out sheepforpca

6. PCA-

PCA的三个文件：

Plink --sheep --file data(生成LD的文件) --extract data (LD).prune.in --recode --out sheepforpca 1sheepforpca.ped 改为5.ped

2sheepforpca.map 改为5.pedsnp

3将sheepforpca 制作成二进制文件输出5b

plink --file hapmap1 --make-bed --out hapmap1

结果为5b.farm即为ped文件的前6列，将5b.farm 改名为5.pedind

Note: 5.pedind 文件中要将第六列-9换成familyID.

参数文件

Genotypename: 5.ped

Snp name: 5.pedsnp

Indivame: 5.pedind

Evecoutname: 5.pca.evec

Evaloutname: 5.eval

Altnormstyle: NO

Numoutevec: 3

Numoutlieriter: 5

Numoutlierevec: 10

Outlier sigmathresh: 6.0

Qt mode: NO

将上述文件拷贝到eigensoft/bin 文件夹内

打开命令

Cd EIG5.01/bin

作图命令

./smartpca –p 5.par

./ploteig –I 5.pca.evec –c 1:2 –p AL:BSB:…Tiberan –x 5-0

即可得到PCA

在5.pca.evec文件中可以看到主成分占的比例。

7 原始SNP数据转化成map和ped文件

>data=read.csv("E:/SNP/zang.csv") (data=read.csv("E:/SNP/chicken.csv")

>Ta=t(data)

> write.table(Ta,file="E:/snp/1.txt", quote=FALSE, sep=" ", na="0")

检查命令：--Compound –genotypes

8: 近交系数计算, 多态性含量, Ho He 哈温P值）

plink --file filename --het --out filename1

plink --file filename --homozyg --out filename1

plink --file filename --hardy --out filename1（Plink --file filename –hardy，结果为plink.hwe）9: ADZE软件计算Ar 和pAr

Plink 转化成structure

1Plink --file filename --recode-structure --out filename1

2用PGDspider 转化ped 文件为structure结构。

将plink转化的位点信息粘贴到PGDspider转化的文件中

全基因组关联分析

plink --file data --remove mylist.txt --recode --out filename

plink --sheep --file filename --out name –assoc

plink --sheep --file filename --out name --assoc --adjust

R软件中绘制manhattan 图

先安装qqman软件

> library(qqman)

> results<- read.table("D:/plink/plink1/plink.assoc",T)

> manhattan(results, ylim = c(0, 10), col = c("blue4", "orange3"))