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(\>Ta=t(data)

> write.table(Ta,file=\检查命令：--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，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(\

> manhattan(results, ylim = c(0, 10), col = c(\

结果为plink.hwe） R软件中绘制qq-plot图 > library(qqman)

> results<- read.table(\> qq(results$P)

用其他关联分析方法：

plink --sheep --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

本文来源：https://www.bwwdw.com/article/hv5x.html