全基因组关联分析（GWAS）：LD Block

画LD block用的比较多的有两个软件，一个是haploview，另一个是R包LDheatmap。

haploview

数据导入

传说plink可以直接转化出HV格式（参数 --recode HV），但是用的时候提示没有case/control表型，不知如何解决。因此自己将plink格式转换为Linkage formate格式
Linkage formate格式有两个文件：ped文件和info文件

ped文件

与plink的ped非常相似：
第一列：家系。没有家系可以指定其它唯一ID
第二列：个体。
第三列、第四列：父本、母本。0 表示缺失
第五列：性别，1=雄性，2=雌性。注意，必须选择一个，没有缺失值
第六列：个体的患病状态，0=未知，1=未患病，2=患病。

info文件

第一列：SNP ID
第二列：SNP 物理位置

转换

plink1.7

# plink --bfile myfile --snp snpID --window 1000 --recode tab --out test
plink --file myfile --from-bp updream --to-bp downdreanm --tab --recode --out test

1. 对于植物来说需要重新设定ped文件第五列性别（如都设为 1）和第六列（都设为0），之后此文件直接作为Linkage Formate的ped文件。
2. 提取map文件的第二列（SNP ID）和第四列（physical）作为info文件

读入数据

可以用图形界面，也可以用命令行

java -jar ~/biotools/Haploview.jar -memory 3000 -pedfile test.ped -info tes.info -skipcheck -svg

若加上-nogui参数，可以不打开图形界面，直接输出LD block图片。
注意：

1. 用图形界面的好处是可以不断删除低质量的SNP（如MAF太小），调整LD block
2. 但是若SNP数量过多，会很耗内存。

LDheatmap

输入文件

• 第一个参数有两种类型，LD matrix（计算好的相关系数矩阵）或者基因型文件（让 LDheatmap 自己计算）
• 第二个参数是代表SNP物理位置的向量。

# 提取目的范围的SNP
plink --file myfile --from-bp updream --to-bp downdreanm --tab --recode --out outfile
# 提取第二个文件所需数据
cut -f4 outfile.map >dist.txt
# 提取SNP ID备用
cut -f2 outfile.map >name.txt
# 用plink计算LD matrix，输出文件为plink.ld
plink -file outfile --r2 --matrix
# 提取基因型文件
cut -f7- outfile.ped >geno.txt
#Alleles 之间`/`分隔
sed -i 's/ /\//g' geno.txt

library(LDheatmap)
library(genetics)
#导入dist.txt和name.txt
dist <- read.table('dist.txt')
name <_ read.table('name.txt')
#第一种方法，导入LD matrix，读入并转换为maxtrix格式
mx <- as.matrix(read.table('plink.ld'))
#矩阵行列名用SNP ID替换，以便作图时显示SNP名字
colnames(mx) <- name[,1]
rownames(mx) <- name[,1]
#作图
ld <- LDheatmap(mx,genetic.distances = dist[,1])
#第二种方法，导入基因型文件
geno <- read.table('geno.txt')
colnames(geno) <- name[,1}
#需要转换成SnpStats的genotype格式
num <- ncol(geno)
for (i in 1:num){geno[,i] <- as.genotype(geno[,i])}
#作图
ld <- LDheatmap(geno, genetic.distances = dist[,1])
#输出结果 ld 可以直接作为输入参数
#用colorRampPalette进行调色 
color = c(colorRampPalette(colors = c("red","white"))(80),colorRampPalette(colors = c("white","grey"))(20))
color2 <- colorRampPalette(c("red","#FA9FB5","#BDBDBD"))
LDheatmap(ld,color = color2(100),SNP.name = c('SNP1','SNP2'))