E:\ > cd e:
E:
E:\ > cd plink-1
E:\plink-1>plink –file test
- Map 更新
Plink --goat --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名 - Quality control
Call rate >98%/99%
Plink --file goat --geno 0.02 --recode --out goatgeno
Plink --file goatgeno --mind 0.01 --recode --out goatmind
MAF>0.05
Plink --file goatmind --maf 0.05 --recode --out goatmaf
Hardy-Weinberg equilibrium <0.0001
Plink --file goatmaf --hwe 0.0001 --recode --out goathwe
Exclude the SNP markers with either chromosome or both unknown
Plink --goat --file goathwe --extract 4newsnp.txt --recode --out goat4
Note: 制作4newsnp.txt（包含 chromosome 和 base-pair position 都为0的SNP）
To identify sample duplication or half-sibs or closer
Plink –goat –file goat4 –genome –max 0.85
Note：Check the genome file - LD quality control
Plink –goat --file goat4 –indep-pairwise 100 25 0.2 –out goatld0.2
Plink --goat --file goat4 --indep-pairwise 100 25 0.05 --out goatld0.05
Plink --file goat4 --ld-window-r2 0.2 --out goatldr0.2
输出结果为 data prunein 和 data prune out
(质控时，要去除X染色体)
将data prune in 转化为ped和map
Plink --goat --file 114hwe --extract 114goat0.05.prune.in --recode --out goatforpca - PCA-
PCA的三个文件：
Plink --goat --file data(生成LD的文件) --extract data (LD).prune.in --recode --out goatforpca
1goatforpca.ped 改为5.ped
2goatforpca.map 改为5.pedsnp
3将goatforpca 制作成二进制文件 输出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 --goat --file filename --out name –assoc
plink --goat --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"))
R软件中绘制qq-plot图
library(qqman)
results<- read.table("D:/plink/plink1/plink.assoc",T)
qq(results$P)
用其他关联分析方法：
plink --goat --file mar --out mar-model --model --model-trend --adjust
LD 分析（haploview）
1 info 文件生成：plink --file hu-M2 --recodeHV --out hu-M2HV
R 安装 GenABEL
Install packages (GenABEL)
Install packages (MASS)
Install packages(Gen ABEL.data)
加载安装包
Library (MASS)
Library(Gen ABEL.data)
Library(GenABEL)
在使用GenABEL前需要准备4个文件
Ped、map、phen(当ped中含有多个表型时用到)、praw
1生成tped、tfram 文件
Plink –file name –transpose –recode –out gwa-gabel
当多个表型时还还需要—pheo phenol.phen –pheno-name
2制作praw文件
格式 id sex phen(sex:female=0, male=1 phen case=1 control=0)
“Sss12” 1 0
“Sss18” 1
1 GLM test
Testb<- scan.glm(‘phen~CRSNP’, family=binomial(),data=b.dat)
2 score test
Testb.qt<- qtscore(phen, data=b.dat, trait=”binomial”)
Test.qt@lambda
格式转化：
1 二进制格式转为 tped/tfam 或 ped/map
用bfile来读取test3.bed，test3.bim和test3.fam文件
生成test4.tped和test4.tfam
plink --bfile test3 --recode --transpose --out test4
生成test5.ped和test5.map
plink --bfile test3 --recode --out test5
- Quality control
Call rate >98%/99%
Plink --file goat1 --geno 0.02 --recode --out goatgeno
Plink --file goatgeno --mind 0.01 --recode --out goatmind
MAF>0.05
Plink --file goatmind --maf 0.05 --recode --out goatmaf
Hardy-Weinberg equilibrium <0.0001
Plink --file goatmaf --hwe 0.0001 --recode --out goathwe
Exclude the SNP markers with either chromosome or both unknown
Plink --file goathwe --exclude 4newsnp.txt --recode --out goat4
Note: 制作4newsnp.txt（包含 chromosome 和 base-pair position 都为0的SNP）
To identify sample duplication or half-sibs or closer
Plink –goat –file goat4 –genome –max 0.85
Note：Check the genome file - LD quality control
Plink --file goat --indep-pairwise 100 25 0.2 --out goatld0.2
Plink --goat --file goat --indep-pairwise 100 25 0.05 --out goatld0.05
Plink --file goat4 --ld-window-r2 0.2 --out goatldr0.2
输出结果为 data prunein 和 data prune out
(质控时，要去除X染色体)
将data prune in 转化为ped和map
Plink --goat --file goat --extract goatld0.05.prune.in --recode --out goatforpca
D:\科研\山羊SNP\SNP分析软件\plink-1.07-dos
plink --cow --vcf Europeanibex.vcf --recode --out output
山羊分群：
