画LD block用的比较多的有两个软件,一个是haploview,另一个是R包LDheatmap。
haploview
数据导入
plink 1.9 可以直接转化出haploview格式(参数 --recode HV)
plink --bfile ../snp --recode HV tab --out snp
每一个染色体都会产生两个文件(Linkage format 和 map info),如, snp.chr-1.ped和snp.chr-1.info,使用方式如下
java -jar Haploview.jar -pedfile snp.chr-1.ped -info snp.chr-1.info
也可以自己将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
- 对于植物来说需要重新设定ped文件第五列性别(如都设为 1)和第六列(都设为0),之后此文件直接作为Linkage Formate的ped文件。
- 提取map文件的第二列(SNP ID)和第四列(physical)作为info文件
读入数据
可以用图形界面,也可以用命令行
java -jar ~/biotools/Haploview.jar -memory 3000 -pedfile test.ped -info tes.info -skipcheck -svg
若加上-nogui
参数,可以不打开图形界面,直接输出LD block图片。
注意:
- 用图形界面的好处是可以不断删除低质量的SNP(如MAF太小),调整LD block
- 但是若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'))