2020.11.9【GWAS/WGS流程】丨全基因組關(guān)聯(lián)分析繪圖全流程_穆易青的博客-CSDN博客_gwas全基因組關(guān)聯(lián)分析

一、數(shù)據(jù)篩選及格式轉(zhuǎn)換

Ｅ跎肌Ｋ鸺蟆恒水！注意：如果要同時使用蜕提，正確做法是先SNP（geno）后樣本(mind)

1.去除多等位基因，indel

nohup /home/software/bcftools/bcftools view -m 2 -M 2 --type "snps"  75.vcf.gz -O z -o 75.recordsnps2vcf.gz &

2.格式轉(zhuǎn)換（map ped to bam bed fam）

nohup/home/software/plink --allow-extra-chr --chr-set 29 --file Neogen_China_CHN_BOVG100V1_20201104 --make-bed --out Neogen_China_CHN_BOVG100V1_20201104 &

3.LD過濾 #（做進化樹需要用到這一步）

# 基于連鎖不平衡的SNP修剪（窗寬50凡傅、刪除LD大于0.2的SNP對中一個速警、將窗口向前移動25個SNP）
nohup /home/software/plink --allow-extra-chr  --chr-set 29 -bfile Neogen_China_CHN_BOVG100V1_20201104  --indep-pairwise  50 25 0.2 --out ld.Neogen_China_CHN_BOVG100V1_20201104-502502 &
# 提取過濾好的
nohup /home/software/plink --allow-extra-chr --chr-set 29 -bfile Neogen_China_CHN_BOVG100V1_20201104 --extract ld.Neogen_China_CHN_BOVG100V1_20201104-502502.prune.in --make-bed --out ld.Neogen_China_CHN_BOVG100V1_20201104-502502 &

4.最終QC進行質(zhì)控，位點過濾（做群體結(jié)構(gòu)需要這一步 ）剔除高缺失率（--geno ）和極低等位基因頻率( --maf )的SNP

nohup /home/software/plink --allow-extra-chr --chr-set 29 -bfile ld.Neogen_China_CHN_BOVG100V1_20201104-502502 --geno 0.05 --maf 0.03 --make-bed --out ld.QC.Neogen_China_CHN_BOVG100V1_20201104-502502-geno005-maf03 &
# --geno 0.05   大于95%的個體都具有的變異位點才保留强霎，其他去除忿项；或者說保留檢出率高于0.95的SNP。
# --maf 0.03    次等位基因頻率脆栋，頻率較低的第二等位基因的頻率（防止假陽性）倦卖；將maf < 0.03的SNP篩選出來并過濾掉，僅包括MAF >= 0.03的SNP椿争，默認(rèn)值為0.03怕膛。
# --mind 0.10   去除基因型丟失率大于10%的個體樣本；
# --hwe 0.0001  保留符合Hardy-Weinbery 的變異位點秦踪。

5.格式轉(zhuǎn)換

# bed轉(zhuǎn)成map ped(待定)
nohup /home/software/plink --allow-extra-chr --chr-set 29  --bfile ld.QC.Neogen_China_CHN_BOVG100V1_20201106-502502-geno02-maf03 --recode --out ld.QC.Neogen_China_CHN_BOVG100V1_20201106-502502-geno02-maf03 &
# 格式轉(zhuǎn)換 將bed bim fam 轉(zhuǎn)vcf文件
nohup /home/software/plink  --allow-extra-chr --chr-set 29 --bfile test --recode vcf-iid --out test &

二褐捻、群體結(jié)構(gòu)之構(gòu)建進化樹（22.10.25 注意數(shù)據(jù)不要ld過濾）

1.iqtree 構(gòu)建進化樹(基于每一個snp畫，較準(zhǔn)確但用時長）

# 轉(zhuǎn)成fa 
nohup python /home/hmy/ped2fa.py   Neogen_China_CHN_BOVG100V1_20201106-geno02-maf03.ped     ld.QC.Neogen_China_CHN_BOVG100V1_20201106-502502-geno02-maf03.fa &
# 過濾I成N 
nohup sed -i "s/I/N/g" Neogen_China_CHN_BOVG100V1_20201104-geno02-maf03.fa &
# mafft比對 
nohup  mafft --auto Neogen_China_CHN_BOVG100V1_20201104-geno02-maf03.fa > Neogen_China_CHN_BOVG100V1_20201104-geno02-maf03_aligned.fa &
# iqtree 計算進化樹
nohup /home/software/iqtree-1.6.12-Linux/bin/iqtree -s Neogen_China_CHN_BOVG100V1_20201104-geno02-maf03_aligned.fa  -m TEST   -st DNA   -bb 1000 -nt AUTO &
# 9.iTOL畫樹( 有時畫圖需替換ID名稱）
https://itol.embl.de/upload.cgi
輸入XXX.bionj 文件

# ped2fa.py腳本
#!/usr/bin/env python
import sys
input_ped = sys.argv[1]
output_fa = sys.argv[2]
a1=open(input_ped , "r")
x = a1.readline()
b1 = open(output_fa , "w")
i=0
while x:
    x = x.strip().split()
    c= (">",x[0])
    columns= ''.join(c)
    b1.write(columns + "\n")
    seq = ''.join(x[6:])
    sequence = seq.replace("0", "N")
    leng =int(len(sequence))
    len1 =leng-1
    all = ("A","T","C","G")
    while i<len1:
        b1.write(sequence[i])
        i = i+2
    b1.write("\n")
    x = a1.readline()
    i=0
a1.close
b1.close

2.孫老師腳本畫（基于親緣關(guān)系畫the neighbor-joining phylogenetic tree was constructed by genetic distance (1-IBS)）

# 畫進化樹ibs
nohup /home/software/plink --allow-extra-chr --chr-set 29  --noweb --file com-328-100k-need_noneed.ped-geno02-maf03 --genome --out com-328-100k-need_noneed.ped-geno02-maf03 &
(使用質(zhì)控后的ped和map文件椅邓，產(chǎn)生com-328-100k-need_noneed.ped-geno02-maf03.genome文件)
# 使用孫老師腳本（注意更改該腳本中的樣品數(shù)Ｄ选！＞澳佟）
nohup perl 3____PLINK_genome_MEGA.pl &
# 使用mega軟件查看進化樹板壮，mega軟件中可刪除樣品看進化樹！：献　４戮！"3-22-lzx_Dis.meg"在mega軟件打開透葛。

#3____PLINK_genome_MEGA.pl腳本：
#!usr/bin/perl
# define array of input and output files

open (AAA,"0-1-cattle-293-chr1-29-snp-indel.filter.pass-geno01-maf005.genome") || die "can't open AAA";
open (BBB,"/home/ysq/20210301-cattle-vcf-hebing/final-vcf/NEW-3-21/1-cattle-293-chr1-29-snp-indel.filter.pass-geno01-maf005.fam") || die "can't open BBB"; 
open (CCC,">3-22-lzx_Dis.meg");
my @aa=<AAA>;
my @bb=<BBB>;

$sample_size=293; ### 樣品個數(shù)

print CCC "#mega\n!Title: $sample_size pigs;\n!Format DataType=Distance DataFormat=UpperRight NTaxa=$sample_size;\n\n"; 

foreach ($num1=0;$num1<=$#bb;$num1++){
    chomp $bb[$num1];
    @arraynum1=split(/\s+/,$bb[$num1]);
    print CCC "#$arraynum1[1]\n";       ##
    }

print CCC "\n";


@array=();
foreach ($num2=1;$num2<=$#aa;$num2++){
    chomp $aa[$num2];
    
    
    @arraynum1=split(/\s+/,$aa[$num2]);
    push(@array,1-$arraynum1[12]);
    }

    
@array2=(0);
$i=$sample_size;
while ($i>0){   
    push(@array2,$array2[$#array2]+$i);
    $i=$i-1;
    }
print "@array2";

for ($i=($sample_size-1); $i>=0; $i=$i-1){
    print CCC " " x ($sample_size-($i+1));
        for ($j=$array2[$sample_size-$i-1]; $j<=$array2[$sample_size-$i]-1; $j++){
                                                                                  print CCC "$array[$j] ";  
                                                                                 }
        print  CCC "\n";
    }

close AAA;
close BBB;
close CCC;

三笨使、群體結(jié)構(gòu)之a(chǎn)dmixture(注意：重測序數(shù)據(jù)需替換染色體NCto1)

！Ａ藕Α硫椰！在構(gòu)建admixture前需要考慮樣本中每個群體中個體的數(shù)目，在總樣本分析好后，與刪除個別個體再次分析后靶草，兩次得到圖形的顏色占比會發(fā)生變化蹄胰。

# 使用bed文件,其中重測序數(shù)據(jù)中SNP沒有名稱，需要補上名稱奕翔，使用awk:
awk  '{print $1,$1":"$4,$3,$4,$5,$6}'  XXX.bim > XXX-change.bim
sed -i 's/\t/ /g'  XXX-change.bim   # 去除vcf文件中的tab鍵
sed -i 's/ /\t/g'  XXX-change.bim   # 空格換成tab鍵
# 畫圖admixture可以使用gz壓縮文件,做出bed
nohup /home/software/plink --allow-extra-chr --chr-set 27 --vcf XXX.vcf.gz  --make-bed --out XXX
# ld過濾,maf過濾等
nohup /home/software/plink --allow-extra-chr  --chr-set 29 -bfile test  --indep-pairwise  50 25 0.2 --out ld.test-502502 &
nohup /home/software/plink --allow-extra-chr --chr-set 29 -bfile test --extract ld.test-502502.prune.in --make-bed --out ld.test-502502 &
#轉(zhuǎn)成map ped 
nohup /home/software/plink --allow-extra-chr --chr-set 27  --bfile test  --recode --out test &
#格式轉(zhuǎn)換 將bed bim fam 轉(zhuǎn)vcf文件 
nohup  /home/software/plink --allow-extra-chr --chr-set 27 --bfile test --export vcf --out test &(.log/ .nosex/ .vcf文件)
#使用bed文件烤送，計算K
##Error：Invalid chromosome code!  Use integers.
# (判斷因為染色體為NC格式，更改成1糠悯、2帮坚、3格式)替換染色體為數(shù)字后，轉(zhuǎn)換為bed文件
# 計算K值
for K in 1 2 3 4 5; do /home/software/admixture_linux-1.3.0/admixture --cv combine_100k-328_hebing-bcftools.bed $K | tee log${K}.out; done
grep -h CV log*.out
#提取vcf文件中的個體list 
nohup /home/software/bcftools/bcftools query -l test.vcf > test-id.txt &
nohup /home/software/vcftools/src/cpp/vcftools --vcf test.vcf --plink --out test & （vcf to map不會影響fid互艾、iid试和，推薦）
nohup /home/software/plink --allow-extra-chr --chr-set 29 --file test --make-bed --out test &（map to bed不會影響fid、iid纫普，推薦）
nohup /home/software/plink  --allow-extra-chr --chr-set 29 --bfile test --recode vcf-iid --out test &（bed to vcf 不會影響fid阅悍、iid，推薦）

制作order.txt文件（也可為.csv文件）
三列： 1.地區(qū)Asia 2.ID名稱 3.樣本品種
1.fam文件中的IID（第二列）更換為自己想要的個體名稱昨稼，可批量替換节视。2.制作一個ld.QC.Neogen_China_CHN_BOVG100V1_20201104-502502-geno02-maf03.order.txt文件(order.txt文件需要另存為制表符分隔格式文件)，可將order.txt文件的第三列改成真正的個體名稱假栓，這樣圖中就會顯示每個個體名稱寻行。
3.可將第一列換成123...
4.表格中無空格
Africa Dp-8 Dp
Africa Dp-9 Dp
asia Han-1 STH
asia Han-10 STH
asia Han-2 STH
asia SRS1227526 Tan
asia SRS1227527 Tan
Middle_East ERS154526 wild
Middle_East ERS154528 wild
Middle_East ERS154529 wild
Middle_East ERS154530 wild

### R腳本畫圖 
##### 1.未排序直接畫圖
tbl=read.table("hapmap3.3.Q")
barplot(t(as.matrix(tbl)), col=rainbow(3),xlab="Individual #", ylab="Ancestry", border=NA)

##### 2.排序后畫圖
#將Q文件，bed文件匾荆，fam文件拌蜘，order.txt文件均放入同一個文件夾下
#!/usr/bin/env Rscript
## sort
sort.admixture <- function(admix.data){
## sort columns according to the cor
    
    k <- length(admix.data)
    n.ind <- nrow(admix.data[[1]])
    name.ind <- rownames(admix.data[[1]])
    admix.sorted <- list()
    
    if (admix.data[[1]][1,1] > admix.data[[1]][1,2]){
        admix.sorted[[1]] <- admix.data[[1]]
    }else{
        admix.sorted[[1]] <- admix.data[[1]][,c(2,1)]
    }
    
    for (i in 1:(k-1)){
        admix <- matrix(nrow = n.ind, ncol = (i + 2))
        cors <- cor(admix.sorted[[i]], admix.data[[i + 1]])
        sorted.loc <- c()
        for (j in 1:nrow(cors)){
            cor <- cors[j,]
            cor[sorted.loc] <- NA
            sorted.loc <- c(sorted.loc, which.max(cor))
        }
        
        sorted.loc <- c(sorted.loc, which(! 1:ncol(cors) %in% sorted.loc))
        cat("n_max = ", sorted.loc, "\n")
        admix <- admix.data[[i + 1]][,sorted.loc]
        rownames(admix) <- name.ind
        admix.sorted[[i + 1]] <- admix
    }
    return(admix.sorted)
}

sort.iid <- function(k.values, groups){
    ##k.values <- admix.values[[1]]
    ##groups <- admix.fam
    max.col <- which.max(colSums(k.values))
    k.values <- cbind(k.values, groups[match(rownames(k.values), as.character(groups$iid)),])
    k.values <- transform(k.values, group = as.factor(k.values$fid))
    k.means <- tapply(k.values[,max.col], k.values$group, mean)
    k.means <- k.means[order(k.means)]
    k.sort <- data.frame(id = names(k.means), 
                         order = order(k.means),
                         mean = k.means)
    k.values$order <- k.sort[match(as.character(k.values$group), k.sort$id), 3]
    k.values <- k.values[order(k.values$order, k.values[,max.col]),]
    return(rownames(k.values))
}

sort.fid <- function(iid.order, fid.order, fam.table){
    new.order <- c()
    for (fid in fid.order){
        new.order <- c(new.order, which(iid.order %in% fam.table[fam.table$fid == fid, "iid"]))
    }
    return(iid.order[new.order])
}

read.structure <- function(file, type = "structure"){
    if (type == "structure"){
        k.values <- read.table(file = file, header = F)
        rownames(k.values) <- k.values[,1]
        k.values[,1:3] <- NULL
    }else{
        k.values <- read.table(file = file, header = F)
    }
    return(k.values)
}

add.black.line <- function(data, groups, nline = 1){
    # data <- as.matrix(plot.data[[2]])
    # groups <- group.name
    # nline <- 3
    
    data <- as.matrix(data)
    group.name <- unique(groups)
    new.data <- matrix(NA, ncol = ncol(data))
    black.data <- matrix(NA, nrow = nline, ncol = ncol(data))
    new.name <- c(NA)
    for (name in group.name){
        new.data <- rbind(new.data, black.data)
        new.data <- rbind(new.data, data[which(groups == name),])
        new.name <- c(new.name, rep(NA,nline))
        new.name <- c(new.name, rownames(data)[which(groups == name)])
    }
    
    added.data <- new.data[(nline + 2):nrow(new.data),]
    rownames(added.data) <- new.name[(nline + 2):nrow(new.data)]
    return(added.data)
}


header <- "ld.com-328-100k-64-50k-need_noneed-502502-geno02-maf03"       
max.k <- 10   ###文件夾下有多少Q(mào)文件，此處數(shù)值就寫幾

admix.fn <- paste(header, 2:max.k, "Q", sep = ".")
fam.fn <- paste(header, "fam", sep = ".")
admix.fam <- read.table(fam.fn, stringsAsFactors = F,
                        col.names = c("fid", "iid", "pid", "mid", "sex", "pheno"))
admix.values <- lapply(admix.fn, read.table, header = F, 
                       row.names = as.character(admix.fam$iid))
order.fn <- paste(header, "order.txt", sep = ".")
admix.order <- read.table(order.fn, col.names = c("region", "iid", "fid"), stringsAsFactors = F)
id.order <- admix.order$iid
#id.order <- sort.iid(admix.values[[1]], admix.order)
admix.data <- list()
for (i in 1:length(admix.values)){
    admix.data[[i]] <- admix.values[[i]][id.order,]
}

species <- as.character(admix.order[,1])

sorted.data <- sort.admixture(admix.data)

## add black line in plot
nline <- 1
plot.data <- list()
group.order <- admix.order[match(id.order, admix.order$iid),3]
for (i in 1:length(sorted.data)){
    plot.data[[i]] <- add.black.line(sorted.data[[i]], group.order, nline = nline)
}

## add xlab to plot
plot.id.list <- rownames(plot.data[[1]])
#plot.xlab <- admix.fam[match(x = plot.id.list, table = admix.fam$iid),]
plot.xlab <- admix.order[match(x = plot.id.list, table = admix.order$iid),]

plot.lab <- unique(plot.xlab$fid)
plot.lab <- plot.lab[!is.na(plot.lab)]
plot.at <- c()
start <- 0
for (fid in plot.lab){
    xlen <- length(which(plot.xlab$fid == fid))
    gap <- start + floor(xlen / 2)
    plot.at <- c(plot.at, gap)
    start <- start + nline + xlen
}


## barplot admixture and structure
library(RColorBrewer)

my.colours <- c(brewer.pal(8, "Dark2"), "mediumblue", "darkred", "coral4",
                "purple3", "lawngreen", "dodgerblue1", "paleturquoise3",
                "navyblue",  "green3", "red1", "cyan", 
                "orange", "blue", "magenta4", "yellowgreen", "darkorange3", 
                "grey60", "black")
max.k <- length(plot.data)
n <- dim(plot.data[[1]])[1]
#pdf(file=paste(header, "admix.plot.pdf", sep = "."), width = 16, height = 12)
png(file=paste(header, "admix.plot.png", sep = "."), res=300, width = 2400, height = 2000)
par(mfrow = c(max.k,1),  mar=c(0.5,2,0,0), oma=c(6,0,1,0))
par(las=2)
#for (i in 1:(max.k - 1)){

for (i in 1:max.k){
    barplot(t(as.matrix(plot.data[[i]])), names.arg = rep(c(""), n), 
            col = my.colours, border = NA, space = 0, axes = F, 
            ylab = "")
    axis(side = 2, at = 0.5, labels = as.character(i+1), tick = F, hadj = 0)
}


axis(side = 1, at = plot.at, labels = plot.lab, tick = F, lty = 15)
#barplot(t(as.matrix(plot.data[[max.k]])), names.arg = rownames(plot.data[[1]]), axes = F,
#        col = my.colours, border = NA, space = 0, 
#        ylab = paste("K=", max.k + 1, sep = ""), cex.axis=0.6, cex.names=0.6)
dev.off()

>畫一個Q文件的腳本
tbl=read.table("/USER/zhusitao/Project/xj/reseq/result/11.admixture/QC.7.Q");
pdf("/USER/zhusitao/Project/xj/reseq/result/11.admixture/Q7.pdf")
barplot(t(as.matrix(tbl)),col= rainbow(7),xlab="Individual", ylab="Ancestry",border = NA,space = 0);
dev.off()

四、群體結(jié)構(gòu)之主成分分析PCA

# 使用常染色體的vcf數(shù)據(jù)進行畫圖（染色體為NC格式，更改成1座硕、2、3格式）
# 生成用于PCA分析的matrix
nohup /home/software/plink --allow-extra-chr --chr-set 27 --vcf ld.test3-502502-chr-33.vcf --make-bed --out ld.test3-502502-chr-33 &
nohup /home/software/gcta_1.92.3beta3/gcta64  --bfile ld.test3-502502-chr-33 --make-grm --autosome-num 27 --out ld.test3-502502-chr-33.gcta &
 # --autosome是只分析常染色體
nohup /home/software/gcta_1.92.3beta3/gcta64  --grm ld.test3-502502-chr-33.gcta  --pca 3 --out ld.test3-502502-chr-33.gcta.out &
# 得到的out.eigenvec即可用于下一步R作圖使用举娩。

制作ld.test3-502502-chr-33.txt表格,共四列,第1列為品種名,24列為樣品名,3列為順序數(shù)字
SAMM MLN-1 1 MLN-1
SAMM MLN-2 3 MLN-2
SAMM MLN-3 4 MLN-3
SAMM MLN-4 5 MLN-4
CM SRR11657659 11 SRR11657659
CM SRR11657660 12 SRR11657660

### R腳本畫圖
 
a=read.table("test.eigenvec",header=F)
head(a)
dim(a)
b=read.table("ld.test3-502502-chr-33.txt",header=F)
head(b)

library("ggplot2")  
qplot(a[,3],a[,5],col=b[,1])
# 比較a文件3和5列,需要比較34,45
Breed=b[,1]
ggplot(data  = a ,
       aes(x = a[,3],
           y = a[,4],   #此處也做相應(yīng)更改
           group = Breed,
           shape = Breed,
           color = Breed)
)+geom_point(size=2) +scale_shape_manual(values = seq(0,75))
# 通常而言，PCA圖會和系統(tǒng)發(fā)育樹以及群體結(jié)構(gòu)一起解釋构罗，相互驗證铜涉。
###PCA.explain.R
m <- as.matrix(read.table("test.eigenval",header=F))
explainm=m/sum(m)
explainm[1:3]  #此處上面步驟 --pca 3，一般3個可以解釋
plot(explain)
PCA方差解釋率顯示
m <- as.matrix(read.table("test.eigenval",header=F))
explainm=m/sum(m)
explainm[1:3]
plot(explain)

五绰播、LD衰減圖

如果要理解LD衰減圖骄噪，我們就必須先理解連鎖不平衡（Linkagedisequilibrium尚困，LD）的概念蠢箩。連鎖不平衡是由兩個名詞構(gòu)成，連鎖+不平衡。從一個類似的概念入手谬泌，更容易理解LD的概念滔韵，那就是基因的共表達。換句話來說掌实，當(dāng)位于某一座位的特定等位基因與另一座位的某一等位基因同時出現(xiàn)的概率大于群體中因隨機分布的兩個等位基因同時出現(xiàn)的概率時陪蜻，就稱這兩個座位處于連鎖不平衡狀態(tài)（linkage disequilibrium）。如果兩個SNP標(biāo)記位置相鄰贱鼻，那么在群體中也會呈現(xiàn)基因型步調(diào)一致的情況宴卖。如果兩個基因座是相關(guān)的，我們將會看到某些基因型往往共同遺傳邻悬，即某些單倍型的頻率會高于期望值症昏。
這種不同基因座間的相關(guān)性，用一個數(shù)值來衡量就是D值父丰。類似相關(guān)系數(shù)是標(biāo)準(zhǔn)化后的協(xié)方差肝谭，LD系數(shù)（r2）則是標(biāo)準(zhǔn)化后的D值，這個數(shù)值在0~1波動蛾扇。r2=0就是兩個位點完全不相關(guān)攘烛，群體中單倍型分布是隨機的（觀測值=期望值）。r2=1就是兩個位點完全相關(guān)镀首，某些基因型（A）只與特定的基因型（B）共同出現(xiàn)坟漱。
一般而言，兩個位點在基因組上離得越近更哄，相關(guān)性就越強靖秩，LD系數(shù)就越大。反之竖瘾，LD系數(shù)越小沟突。也就是說，隨著位點間的距離不斷增加捕传，LD系數(shù)通常情況下會慢慢下降惠拭。這個規(guī)律，通常就會使用LD衰減圖來呈現(xiàn)庸论。LD衰減圖就是利用曲線圖來呈現(xiàn)基因組上分子標(biāo)記間的平均LD系數(shù)隨著標(biāo)記間距離增加而降低的過程职辅。

六、Haploview-連鎖不平衡分析

當(dāng)位于某一座位的特定等位基因與另一座位的某一等位基因同時出現(xiàn)的概率大于群體中因隨機分布的兩個等位基因同時出現(xiàn)的概率時聂示，就稱這兩個座位處于連鎖不平衡的狀態(tài)域携。目前，連鎖不平衡分析是群體進化的經(jīng)典分析條目鱼喉，分析的軟件主要有plink和haploview秀鞭。Haploview是一款分析單倍型的軟件趋观，依托于java形成的可視化界面，以下是整理的相關(guān)介紹及使用方法锋边。

安裝JAVA皱坛、配置環(huán)境，安裝Haploview
提取位置文件
轉(zhuǎn)換成map ped
將vcf文件轉(zhuǎn)換為map ped文件豆巨、并制作info文件
awk  '!/^#/{print $2"\t"$4 }' test.map > test.info
在軟件中輸入ped info文件后剩辟，出現(xiàn)錯誤：invalid affected status
haploview出現(xiàn)"invalid affected status"的解決方法 haploview彈出這種錯誤是因為haploview的缺失值默認(rèn)為0，而plink文件的缺失值一般用”-9“表示往扔，當(dāng)ped文件的缺失值為”-9“時贩猎，haploview不認(rèn)得該字符，就彈出錯誤提示萍膛。 解決方法為將文件中所有的”-9“替換為0融欧。
解決：在linux界面中，用vim命令編輯ped文件后卦羡，輸入以下字符
sed -i  's/-9/0/g' test.ped
”/-9/0/“表示將該文件的”-9“字符替換為”0“字符噪馏，g表示第一行到最后一行所有的”-9“字符都要進行替換。
在軟件中輸入ped與map文件