題目來自生信技能樹論壇

image.png

統(tǒng)計人類參考基因組的每條染色體長度，每條染色體N的含量，GC含量
因為我下載的是hg38,所以這里用hg38的版本

life is short,I use python
所以我還是主要用python來做
再就是用R，這是我第二想掌握的
最后用shell來做

這個模式也是接下來做其他題目的模式

1.python
首先是我自己寫的腳本，這個腳本不太好看，但還是好用的

import sys
import collections
args=sys.argv
filename=args[1]
count_all=collections.OrderedDict() #構建有順序的字典
Length=0
N_count=0
GC_count=0
print("chr","N_count","GC_count","Length",sep="\t")
for line in open(filename):
   if line.startswith(">"):
      if Length:               #N和GC可能不一定統(tǒng)計的到，但是length肯定是有的
         result=str(N_count)+"\t"+str(GC_count)+"\t"+str(Length) #將所有的結果建立成一個以\t分割的字符串，以該字符串和ID構建字典
         count_all[ID]=result
         N_count=0
         GC_count=0
         Length=0
      ID=line[1:].strip()
   else:
      line=line.upper()
      N_count+=line.strip().count("N")
      GC_count+=line.strip().count("G")+line.strip().count("C")
      Length+=len(line.strip())
result=str(N_count)+"\t"+str(GC_count)+"\t"+str(Length)
count_all[ID]=result
for ID,result in count_all.items():
   print(ID,result,sep="\t")

之所以不好看，主要還是我之前不會用一個key對應value，并且這個value中又形成多個鍵值對的用法，比如：

chr_1
 {'c': 7, 'n': 4, 't': 10, 'a': 13, 'g': 10}
chr_2
 {'c': 5, 'n': 4, 't': 6, 'a': 11, 'g': 8}
chr_3
 {'c': 3, 'n': 4, 't': 6, 'a': 10, 'g': 10}
chr_4
 {'c': 2, 'n': 3, 't': 4, 'a': 9, 'g': 7}

所以接下來給出這種用法的代碼
代碼很好看，語法很簡潔，速度也很快，這就是我比較喜歡python的原因

import sys
import collections
args=sys.argv
filename=args[1]
count_ATCG=collections.OrderedDict() #構建有順序的字典
Bases=["a","t","g","c","n"]
for line in open(filename):
   if line.startswith(">"):
      chr_id = line[1:]
      count_ATCG[chr_id] = {}
      for base in Bases:
        count_ATCG[chr_id][base] = 0 ##字典中，chr_id屬于key，base和0共同構成#value，而value中，base又成了key，0成了value
   else:
      line = line.lower()
      for base in Bases:
        count_ATCG[chr_id][base] += line.count(base)

for chr_id, ATCG_count in count_ATCG.items():
    count_sum = sum(ATCG_count.values())
    count_GC = ATCG_count['g'] + ATCG_count['c']
    print(chr_id)
    for base in Bases:
        print("%s: %s" % (base, ATCG_count[base])) #注意占位符的使用

    print("Sum: %s" % count_sum)
    print("N %%: %.2f%%" % (ATCG_count['n']*100.00/count_sum))
    print("GC %%: %.2f%%" % (count_GC*100.00/count_sum))

2.R
先給出R的腳本

seq <- readLines("~/WGS_new/input/genome/new/hg38.fa")  
# 逐行讀取文件

is_ID <- regexpr("^>",seq,perl=T)  
# 對多有文件進行操作，匹配到>開頭的定位1，不是則為-1, 1是指匹配到"^>"的位置，第一個位置

seq_ID <- seq[which(is_ID==1)]
# 將seq中is_ID為1的取出來，也就是將ID行取出來

seq_content <- seq[which(is_ID==-1)]

start <- which(is_ID==1) 
#取出is_ID=1的位置


end <- start[2:length(start)]-1
#取出除了最后一個元素外，其他的每個非ID元素的最后一個位置

end <- c(end,length(seq))
# 把最后一個位置補上

distance <- end - start
# 這算的是每一個ID對應的多少行序列
seq_num_position <- 1:length(start) 
# 算出有多少個ID，按1234...排序

index <- rep(seq_num_position,distance)
# 構建tapply中的因子，用于分組用

seq_content <- tapply(seq_content,index,paste,collapse="")
# 將同一組的序列合并在一起
seq_content <- as.character(seq_content)
# seq_content本來是數(shù)組,as.character()之后轉換成了向量
seq_length <- nchar(seq_content)
# 統(tǒng)計字符串的長度，也就是對應序列的長度


# 計算GC含量和N含量
G_count=""
C_count=""
N_count=""


for ( i in 1:length(seq_content))
{   G_count[i]<-length(gregexpr("[Gg]",seq_content,perl = T)[[i]])+length(gregexpr("[Cc]",seq_content,perl = T)[[i]])
    N_count[i]<-length(gregexpr("[Nn]",seq_content,perl = T)[[i]])
}

#gregexpr與regexpr不一樣，會將所有的匹配的位置輸出


result <- data.frame(seq_ID,G_count,N_count,seq_length)

這里先總結幾個知識點
1.R中逐行讀取數(shù)據(jù)用的是readLines（）
2.which返回的是判斷結果的坐標
3.tapply()函數(shù)，tapply(x,f,g)
需要向量 x (x不可以是數(shù)據(jù)框)，因子或因子列表 f 以及函數(shù) g
tapply（）執(zhí)行的操作是：暫時將x分組，每組對應一個因子水平，得到x的子向量，然后這些子向量應用函數(shù) g
舉個簡單的例子

a <- c(24,25,36,37)
b <- c('q', 'w', 'q','w')
tapply(a, b, mean)
 q  w 
30 31 

4.nchar()函數(shù)用來計算字符串的長度
5.regexpr()函數(shù)
匹配的函數(shù)，返回的是匹配的第一個位置
舉個例子

string <- "abbcbcbcbdbdbdbcbcb" #定義一個字符串
regexpr("c",string,perl=T)    #在string中匹配“c”，返回的是匹配的第一個“c”的位置

# [1] 4
# attr(,"match.length")
# [1] 1
# attr(,"index.type")
# [1] "chars"
# attr(,"useBytes")
# [1] TRUE

如果是匹配一個不存在字符或字符串，返回值為-1

6.gregexpr()函數(shù)
相比regexpr()函數(shù)，gregexpr()函數(shù)表面上就是多了一個“g”,我沒有查，意思應該就是global，也就是把所有的匹配的位置都給返回出來，這樣的話，用于測序的序列，就可以求個length就能算出匹配了多少個，也就是匹配字符的數(shù)量
舉個例子

string <- "abbcbcbcbdbdbdbcbcb"
gregexpr("c",string,perl=T)
# [[1]]
# [1]  4  6  8 16 18
# attr(,"match.length")
# [1] 1 1 1 1 1
# attr(,"index.type")
# [1] "chars"
# attr(,"useBytes")
# [1] TRUE

length(gregexpr("c",string,perl=T)[[1]])
#[1] 5

我不知道有沒有什么其他的方法，但是這個R腳本處理大的數(shù)據(jù)集會很吃力。
原因很簡單這里面很多步驟都重復遍歷了整個文件再做出處理，所以會很慢。而python寫的腳本，是逐行讀取并處理，所以很快，python大概就2~3min就可以處理完，而這個R腳本跑了30多分鐘還沒有結束。

image.png

3.shell(awk)
嘗試用shell主要還是強化一下相關的知識，而且awk我只能算出length
幾個概念:
R:record 行
F:field 列
NR:number of record 行的數(shù)目
NF:number of field 列的數(shù)目
RS:record split 行分割
FS:field split 列分割

time awk 'BEGIN{RS=">";FS="\n"}{if(NR>1){seq="";for(i=2;i<=NF;i++) seq=seq$i; print $1"\t"length(seq)}}' ~/WGS_new/input/genome/new/hg38.fa >count_awk.txt
#seq=""是因為每次循環(huán)的時候，seq值要初始化為""

這個跑完花了2min36s59