R语言中统计基因型数据中每行中各种碱基的数目

 

001、

dir()
dat <- read.table("outcome.ped")
temp <- dat
temp
temp$id1 <- "A"
temp$id2 <- "C"
temp$id3 <- "G"
temp$id4 <- "T"
temp <- t(temp)

result <- matrix(nrow = ncol(temp), ncol = 4)
for (i in 1:ncol(temp)) {
  result[i,] <- rle(sort(temp[,i]))$lengths         ## 统计每列中各种字符的次数
}

result <- as.data.frame(result - 1)
names(result) <- c("A", "C", "G", "T")
dat <- cbind(dat, result)       ## 跟原始数据合并
dat

 

002、pyhton实现

root@PC1:/home/test3# ls
outcome.ped  test.py
root@PC1:/home/test3# cat outcome.ped
G G A A T G
G G G C T T
G G C C T T
G G A C G G
G G C C G G
G G A C T G
root@PC1:/home/test3# cat test.py   ## 实现脚本
#!/usr/bin/python
in_file = open("outcome.ped", "r")
out_file = open("result.txt", "w")
lines = in_file.readlines()
for i in lines:
    temp = i.split()
    A = 0; C = 0; G = 0; T = 0; temp_list = []
    for j in range(len(temp)):
        if temp[j] == "A":
            A = A + 1
        else:
            if temp[j] == "C":
                C = C + 1
            else:
                if temp[j] == "G":
                    G = G + 1
                else:
                    T = T + 1
    temp_list.extend([str(A),str(C),str(G),str(T)])
    out_file.write(' '.join(temp_list) +  "\n")
in_file.close()
out_file.close()
root@PC1:/home/test3# python test.py
root@PC1:/home/test3# ls
outcome.ped  result.txt  test.py
root@PC1:/home/test3# paste outcome.ped result.txt     ## 结果文件
G G A A T G     2 0 3 1
G G G C T T     0 1 3 2
G G C C T T     0 2 2 2
G G A C G G     1 1 4 0
G G C C G G     0 2 4 0
G G A C T G     1 1 3 1