# 生信编程实战第2题（python、R、shell）

image.png

life is short,I use python

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可能不一定统计的到，但是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")

``````

``````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}

``````

``````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

``````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本来是数组,as.character()之后转换成了向量
seq_length <- nchar(seq_content)
# 统计字符串的长度，也就是对应序列的长度

# 计算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)
``````

2.which返回的是判断结果的坐标
3.tapply()函数，tapply(x,f,g)

tapply（）执行的操作是：暂时将x分组，每组对应一个因子水平，得到x的子向量，然后这些子向量应用函数 g

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

4.nchar()函数用来计算字符串的长度
5.regexpr()函数

``````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
``````

6.gregexpr()函数

``````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
``````

image.png

3.shell(awk)

R:record 行
F:field 列
NR:number of record 行的数目
NF:number of field 列的数目
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=""是因为每次循环的时候，seq值要初始化为""
``````