根据位点的物理位置、重组率计算遗传距离

 

001、shell

root@PC1:/home/test2/test# ls
genetic_map.txt
root@PC1:/home/test2/test# head genetic_map.txt | column -t
position  COMBINED_rate(cM/Mb)  Genetic_Map(cM)
61795     0.3516610754          0
63231     0.3500036909          0.0005026053001324
63244     0.3494018702          0.000507147524445
63799     0.3501262382          0.000701467586646
64150     0.3558643956          0.0008263759895016
64934     0.3567249058          0.0011060483156488
65288     0.3633379498          0.001234669949878
66370     10.0482361599         0.0121068614748898
68749     10.051710916          0.0360198817440538
root@PC1:/home/test2/test# awk 'BEGIN{a = 61795; b = 0} NR > 1 {print $0,(($1 - a)/1000000 * $2 + b); a = $1; b = $3}' genetic_map.txt | head | column -t
61795  0.3516610754   0                   0
63231  0.3500036909   0.0005026053001324  0.000502605
63244  0.3494018702   0.000507147524445   0.000507148
63799  0.3501262382   0.000701467586646   0.000701468
64150  0.3558643956   0.0008263759895016  0.000826376
64934  0.3567249058   0.0011060483156488  0.00110605
65288  0.3633379498   0.001234669949878   0.00123467
66370  10.0482361599  0.0121068614748898  0.0121069
68749  10.051710916   0.0360198817440538  0.0360199
69094  10.0518361874  0.0394877652287068  0.0394878

 

 

002、R

dir()
dat <- read.table("testmap.txt", header = F)
dim(dat)
head(dat)
a <- dat[1,1]
b <- 0
result <- vector()
for (i in 1:nrow(dat)) {
  result[i] <- (dat[i,1] - a)/1000000 * dat[i,2] + b
  a <- dat[i,1]
  b <- dat[i,3]
}
dat$cal <- result
head(dat)

 

3、python

root@PC1:/home/test3# ls
genetic_map.txt  test.py
root@PC1:/home/test3# head -n 3 genetic_map.txt
61795 0.3516610754 0
63231 0.3500036909 0.0005026053001324
63244 0.3494018702 0.000507147524445
root@PC1:/home/test3# cat test.py
#!/usr/bin/python

in_file = open("genetic_map.txt", "r")
out_file = open("result.txt", "w")

lines = in_file.readlines()

a = lines[0].split()[0]
b = lines[0].split()[2]
for i in lines:
    out_file.write(i.strip())
    out_file.write("\t")
    i = i.split()
    result = (float(i[0]) - float(a))/1000000 * float(i[1]) + float(b)
    a = float(i[0])
    b = float(i[2])
    out_file.write(str(result))
    out_file.write("\n")

in_file.close()
out_file.close()
root@PC1:/home/test3# python test.py
root@PC1:/home/test3# ls
genetic_map.txt  result.txt  test.py
root@PC1:/home/test3# head -n 3 result.txt
61795 0.3516610754 0    0.0
63231 0.3500036909 0.0005026053001324   0.0005026053001324
63244 0.3494018702 0.000507147524445    0.000507147524445