美文网首页生物信息学与算法Cook R
生信编程实战第2题(python、R、shell)

生信编程实战第2题(python、R、shell)

作者: 天秤座的机器狗 | 来源:发表于2018-08-11 21:57 被阅读115次

    题目来自生信技能树论坛

    image.png

    统计人类参考基因组的每条染色体长度,每条染色体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可能不一定统计的到,但是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本来是数组,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)
    

    这里先总结几个知识点
    1.R中逐行读取数据用的是readLines()
    2.which返回的是判断结果的坐标
    3.tapply()函数,tapply(x,f,g)
    需要向量 x (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
    

    如果是匹配一个不存在字符或字符串,返回值为-1

    6.gregexpr()函数
    相比regexpr()函数,gregexpr()函数表面上就是多了一个“g”,我没有查,意思应该就是global,也就是把所有的匹配的位置都给返回出来,这样的话,用于测序的序列,就可以求个length就能算出匹配了多少个,也就是匹配字符的数量
    举个例子

    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脚本处理大的数据集会很吃力。
    原因很简单这里面很多步骤都重复遍历了整个文件再做出处理,所以会很慢。而python写的脚本,是逐行读取并处理,所以很快,python大概就2~3min就可以处理完,而这个R脚本跑了30多分钟还没有结束。


    image.png

    3.shell(awk)
    尝试用shell主要还是强化一下相关的知识,而且awk我只能算出length
    几个概念:
    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值要初始化为""
    

    这个跑完花了2min36s59

    相关文章

      网友评论

      • 王诗翔:r不会比python慢。你想想用python解决这个问题的,r也可以用同样的思路解决。
        天秤座的机器狗:是的,我R一直不太熟悉
        王诗翔:https://stackoverflow.com/questions/12626637/reading-a-text-file-in-r-line-by-line

      本文标题:生信编程实战第2题(python、R、shell)

      本文链接:https://www.haomeiwen.com/subject/mahdbftx.html