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位點(diǎn)
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
輸出結(jié)果為 data prunein 和 data prune out
(質(zhì)控時(shí)纯陨，要去除X染色體)
將data prune in 轉(zhuǎn)化為ped和map
Plink --goat --file 114hwe --extract 114goat0.05.prune.in --recode --out goatforpca - PCA-
PCA的三個(gè)文件：
Plink --goat --file data(生成LD的文件) --extract data (LD).prune.in --recode --out goatforpca
1goatforpca.ped 改為5.ped
2goatforpca.map 改為5.pedsnp
3將goatforpca 制作成二進(jìn)制文件 輸出5b
plink --file hapmap1 --make-bed --out hapmap1
結(jié)果為5b.farm即為ped文件的前6列偿衰， 將5b.farm 改名為5.pedind
Note: 5.pedind 文件中要將第六列-9換成familyID.
參數(shù)文件
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 文件夾內(nèi)
打開(kāi)命令
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數(shù)據(jù)轉(zhuǎn)化成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: 近交系數(shù)計(jì)算, 多態(tài)性含量, 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女蜈， 結(jié)果為plink.hwe）
9: ADZE軟件計(jì)算Ar 和 pAr
Plink 轉(zhuǎn)化成structure
1Plink --file filename --recode-structure --out filename1
2用PGDspider 轉(zhuǎn)化ped 文件為structure結(jié)構(gòu)虐拓。
將plink轉(zhuǎn)化的位點(diǎn)信息粘貼到PGDspider轉(zhuǎn)化的文件中
全基因組關(guān)聯(lián)分析
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)
用其他關(guān)聯(lián)分析方法：
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前需要準(zhǔn)備4個(gè)文件
Ped沼沈、map、phen(當(dāng)ped中含有多個(gè)表型時(shí)用到)嘱函、praw
1生成tped甘畅、tfram 文件
Plink –file name –transpose –recode –out gwa-gabel
當(dāng)多個(gè)表型時(shí)還還需要—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
格式轉(zhuǎn)化：
1 二進(jìn)制格式轉(zhuǎn)為 tped/tfam 或 ped/map
用bfile來(lái)讀取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
輸出結(jié)果為 data prunein 和 data prune out
(質(zhì)控時(shí)实夹，要去除X染色體)
將data prune in 轉(zhuǎn)化為ped和map
Plink --goat --file goat --extract goatld0.05.prune.in --recode --out goatforpca
D:\科研\(zhòng)山羊SNP\SNP分析軟件\plink-1.07-dos
plink --cow --vcf Europeanibex.vcf --recode --out output
山羊分群：
