搜索
楼主: Jimmy

生信编程直播第二题-hg19基因组序列的一些探究

  [复制链接]

0

主题

5

帖子

115

积分

注册会员

Rank: 2

积分
115
发表于 2017-1-12 20:09:55 | 显示全部楼层
本帖最后由 033-Perl-R-Pyth 于 2017-1-17 18:34 编辑

看视频前:
思路:用测试数据。 制定分割符分割不同的染色体序列,正则 tr 替换
[Perl] 纯文本查看 复制代码
#!usr/bin/perl -w
# Caculate the N, GC content of test hg19
use strict; use warnings;

open A, 'ZY2_test.txt'or die $!;
open OUT, '>Out_ZY2.txt';


$/="\n>";
my ($DNA,$N,$GC, $AT,$len);
$DNA=0;
while (<A>){
        chomp;
        $_=~/(chr_.).*/;
       
        $N=($DNA=~ s/[^ATGC]//ig);
        $GC=($DNA=~ s/[GC]//ig);
        $AT=($DNA=~s/[AT]//ig);
        $len=$AT+$GC;
        
        
        print OUT ">$1\t$len\t$N\t$GC\n";
        print ">$1\t$len\t$N\t$GC\n";
        
}
close A;


运行结果错误。准备看视频学习。
看视频后,修改代码。
[AppleScript] 纯文本查看 复制代码
#!usr/bin/perl -w
# Caculate the N, GC content of test hg19
use strict; use warnings;

open A, 'ZY2_test.txt'or die $!;
open OUT, '>Out_ZY2.txt';


$/="\n>";
my ($chr, $N,$G,$C, $A,$T,$len,$GC);

while (<A>){
	chomp;
	s/>//g;
	($chr)=$_=~/^(chr_.).*/;
	
	$N+=($_=~s/N//ig);
	$G+=($_=~s/G//ig);
	$C+=($_=~s/C//ig);
	$A+=($_=~s/A//ig);
	$T+=($_=~s/T//ig);
	$GC=($G+$C)/($A+$T+$G+$C);
	$len=$A+$T+$G+$C;
		
	print OUT ">$chr\t$len\t$N\t$GC\n";
	print ">$chr\t$len\t$N\t$GC\n";
	
}
close A;


运行结果正确


回复 支持 反对

使用道具 举报

0

主题

7

帖子

136

积分

注册会员

Rank: 2

积分
136
发表于 2017-1-12 20:38:29 | 显示全部楼层
naturehunger 发表于 2017-1-12 19:27
你用的是python2.x还是python3.x版本?这两个版本的除法有一点不同。具体看以下说明和试验。

我的是Python 2.7,谢谢。
回复 支持 反对

使用道具 举报

6

主题

23

帖子

201

积分

中级会员

Rank: 3Rank: 3

积分
201
发表于 2017-1-12 20:55:32 | 显示全部楼层
本帖最后由 qin_qinyang 于 2017-1-12 21:02 编辑

123_python_qin
[AppleScript] 纯文本查看 复制代码
from collections import OrderedDict
seq = ''
chr_d = OrderedDict()
with open("test.fa","r") as f:
    for line in f:
        #line = line.rstrip()
        line=line.strip()
        if line.startswith(">"):
            name = line
        else:
            seq += line.upper()
            seqLen = len(seq)
            N = seq.count("N")
            GC = seq.count("G") + seq.count("C") 
            print (name)
            print (N)
    

>chr_1
3
>chr_1
4
>chr_2
7
>chr_2
8
>chr_3
11
>chr_3
12
>chr_4
15
In [ ]:
结果就是每行都打印N,并叠加到下一行,所以结果并不是每行的N的数目
按照视频修改后,
[Python] 纯文本查看 复制代码
from collections import OrderedDict
seq = ''
chr_d = OrderedDict()
with open("test.fa","r") as f:
    for line in f:
        #line = line.rstrip()
        line=line.strip()
        if line.startswith(">"):
            chr_d[seq] = ''
        else:
            line = line.upper()
            chr_d[seq] += line
    for name,chr_seq in chr_d.items(): #此方法返回元组对的列表
        #print (name,chr_seq)
        seqLen = len(chr_seq)
        N = chr_seq.count("N")
        GC = chr_seq.count("G") + chr_seq.count("C") 
        print (name)
        print (seqLen,N/seqLen,GC/(seqLen-N))

结果如下,就正确了
>chr_1
44 0.09090909090909091 0.425
>chr_2
34 0.11764705882352941 0.43333333333333335
>chr_3
33 0.12121212121212122 0.4482758620689655
>chr_4
25 0.12 0.4090909090909091

回复 支持 反对

使用道具 举报

0

主题

2

帖子

32

积分

新手上路

Rank: 1

积分
32
发表于 2017-1-12 22:08:49 | 显示全部楼层
175 perl

[Perl] 纯文本查看 复制代码
#!/usr/bin/perl
while(<>){
    if(/^>/){
        $chr = $_;
    }
    else{
        $count_A{$chr} += ($_ =~ tr/Aa//);
        $count_T{$chr} += ($_ =~ tr/Tt//);
        $count_C{$chr} += ($_ =~ tr/Cc//);
        $count_G{$chr} += ($_ =~ tr/Gg//);
        $count_N{$chr} += ($_ =~ tr/Nn//);
    }
}
foreach (sort keys %count_A){
    $all_len = $count_G{$_} + $count_C{$_} + $count_A{$_} + $count_T{$_} + $count_N{$_};
    $GC = ($count_C{$_} + $count_G{$_})/$all_len;
    print "$_\t$count_A{$_}\t$count_T{$_}\t$count_C{$_}\t$count_G{$_}\t$GC\n";
}


习惯用s///ig进行替换,很少用tr,去查了下,还是tr最方便。区别如下:
    s///: 替换运算符。s/searchpattern/repalcement/;默认搜索$ _,找出searchpattern,并且用replacement来替换整个匹配的正则表达式。该运算符返回匹配的数量或进行替换的数量,如果没有进行 任何匹配,则返回0。
当然也可以指定某个变量,如$perl=~s/searchpattern/repalcement/; 另外,替换运算符也可以使用非斜杠(/)的界限符,如s###.
    tr/// :转换操作符。tr///的作用与替换运算符有些类似,不过它并不使用正则表达式,转换操作符的句法如下所示:tr /searchment/repalcement/;是把searchment的第一个字符换成replacement的第一个字 符,searchment的第二个字符换成replacement的第二个字符,类推。
回复 支持 反对

使用道具 举报

1

主题

8

帖子

188

积分

版主

Rank: 7Rank: 7Rank: 7

积分
188
发表于 2017-1-12 22:52:19 | 显示全部楼层
写了一个脚本,输入fasta文件,输出'count_atgc.txt'用来计算A,T,G,C的含量和序列长度
[Python] 纯文本查看 复制代码
#!/usr/bin/env python3
#
# Usage example:
#     count_atgc.py tmp.fasta
# 
# Then, you can open the file "atgc_count.txt"

import sys
args = sys.argv

from collections import OrderedDict

chr_dict = OrderedDict()
temp_chr = '' 

with open(args[1], 'r') as f:
	for line in f:
		line = line.strip()
		if line.startswith('>'):
			temp_chr = line
			chr_dict[temp_chr] = ''
		else:
			chr_dict[temp_chr] += line

with open('atgc_count.txt', 'at') as outcome:
    print('seqID\tseqLength\tcount_A\tcount_T\tGC_content', file=outcome)
    for seqName, seq in chr_dict.items():
        seqLen = len(seq)
        seqID = seqName[1:]
        A = seq.count('A') + seq.count('a')
        T = seq.count('T') + seq.count('t')
        G = seq.count('G') + seq.count('g')
        C = seq.count('C') + seq.count('c')
        N = seq.count('N')
        GC = seq.count('G') + seq.count('g') + seq.count('C') + seq.count('c')
        print('{0}\t{1}\t{2}\t{3}\t{4}'.format(seqID, seqLen, A, T, '%.3f'%(GC/(seqLen-N))), file=outcome)


得到的文本结果如下图


本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 1 反对 0

使用道具 举报

0

主题

5

帖子

50

积分

注册会员

Rank: 2

积分
50
发表于 2017-1-12 23:03:45 | 显示全部楼层
本帖最后由 lilith098 于 2017-1-15 19:22 编辑

python+137
状态:还没有看视频,根据基础内容和第一题学到的内容尝试解答。

解题思路:
1 下载hg19的数据,因计算机配置原因,只选择chrY.fa 进行统计;
2 删除序列的空格
3 忽略开头字符“>”的行;
4 将所有列输出到一个字符串中;
5 针对字符串完成计数和输出。

[Python] 纯文本查看 复制代码
import re
import os
from collections import OrderedDict
from operator import itemgetter

os.chdir('D:/Programing/study/biotrainee/test2')

seq = str()
with open('chrY.fa',"r") as chrY:
    for line in chrY:
        line = line.strip()     #去除空格
        if line.startswith('>'):
            chr = line
        else:
            seq += line    #将每一行的字符串相加
#seqLen = len(seq)
N = seq.count('N')
GC = seq.count('G') + seq.count('g') + seq.count('C') + seq.count('c')
print (chr, seqLen, N, GC)

输出报错:
TypeError                                 Traceback (most recent call last)<ipython-input-34-a1cd4a2cd5f8> in <module>()----> 1 seqLen = len(seq)      2 N = seq.count('N')      3 GC = seq.count('G') + seq.count('g') + seq.count('C') + seq.count('c')TypeError: 'int' object is not callable
分步输出:



结果:seqLen = len(seq) 此步骤报错,报错内容:TypeError: 'int' object is not callable

求教各位,这个地方报错的原因是什么?我实在没有看出来。

另外,本次作业只针对一个chrY.fa,没有考虑到题目中test.fa的类型。下一步将尝试处理此问题。

20170115
学习了超长的视频,代码如下:
[AppleScript] 纯文本查看 复制代码
#!/usr/bin/python
#Filename:test_1.py

import os
os.chdir('D:/Programing/study/biotrainee/test2')
count_ATCG = {}
Bases = ['a', 't', 'c', 'g', 'n']

#count the base and store in dict. the key is chr, the value is another dict stored the base number respectively
# with open(sys.argv[1], 'r') as fasta:

with open("test.fa", "r") as fasta:
    for line in fasta:
        if line.startswith('>'):
            chr_id = line[1:]
            count_ATCG[chr_id] = {}
            for base in Bases:
                count_ATCG[chr_id][base] = 0
        else:
            line = line.lower()
            for base in Bases:
                count_ATCG[chr_id][base] += line.count(base)

#print the dict
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))


pycharm运行结果:

chr_1

a: 13
t: 10
c: 7
g: 10
n: 4
Sum: 44
N %: 9.09%
GC %: 38.64%
chr_4

a: 9
t: 4
c: 2
g: 7
n: 3
Sum: 25
N %: 12.00%
GC %: 36.00%
chr_2

a: 11
t: 6
c: 5
g: 8
n: 4
Sum: 34
N %: 11.76%
GC %: 38.24%
chr_3

a: 10
t: 6
c: 3
g: 10
n: 4
Sum: 33
N %: 12.12%
GC %: 39.39%

小结:
1 sys.argv[]输入方式
2 循环语句巩固
3 字典的应用
4 2位浮点百分数的输出

问题:
为什么我的输出不是按照正常的顺序输出?




本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 1 反对 0

使用道具 举报

0

主题

6

帖子

46

积分

新手上路

Rank: 1

积分
46
发表于 2017-1-12 23:27:16 | 显示全部楼层
Python-60

hg19的fa下载了两天都没有下载成功,所以先用测试数据进行了运行,结果如下:

chr_1 44 0.09 0.39
chr_2 34 0.12 0.38
chr_3 33 0.12 0.39
chr_4 25 0.12 0.36

代码复制如下:

[Python] 纯文本查看 复制代码
import collections
import re

path = '/Users/sunshine/practice/practice2-hg19进组序列的一些探究/'

f = open(path+'hg19_test.csv','r')
lines = f.readlines()

chr_dict = collections.OrderedDict()

temp = ''

for line in lines:
    line = re.sub('\>|\n','',line)
    if line.startswith('chr'):
        temp = line
        chr_dict[temp] = ''
    else:
        chr_dict[temp] += line ##纳闷为什么在这里不用定义temp就可以将对应的碱基信息匹配到对应的染色体中,后来想了想,是因为在读取数据的时候,是先读取染色体编号信息的,所以,直接就有染色体的信息在字典中,所以不用重新去匹配,好用!

for chr_name, strand in chr_dict.items():
    length = len(strand)
    N = round((strand.count('N') + strand.count('n'))/length,2)
    GC = round((strand.count('G') + strand.count('C') + strand.count('g') + strand.count('c'))/length,2)
    print (chr_name,len(strand),N,GC)
    
回复 支持 反对

使用道具 举报

0

主题

6

帖子

125

积分

注册会员

Rank: 2

积分
125
发表于 2017-1-13 10:28:30 | 显示全部楼层
023_R
第二题我渣电脑只跑了demo数据,但是收获还是很大的。我看R关于第二题没人贴代码,我贴一下好了
[AppleScript] 纯文本查看 复制代码
seq<-readLines("demo.txt")
is.anno<-regexpr("^>",seq,perl = T)#载入数据保存为向量,并按照“>”标记
is.anno
# [1]  1 -1 -1  1 -1 -1  1 -1 -1  1 -1
mychr<-seq[which(is.anno==1)]
myseq<-seq[which(is.anno==-1)]
start<-which(is.anno==1)
end<-start[2:length(start)]-1
end<-c(end,length(seq))#定位起始终止位置
distance<-end-start
index<-1:length(start)
index<-rep(index,distance)
myseq<-tapply(myseq,index,paste,collapse="")#按照分组整合
seq.content<-as.character(myseq)
seq.len<-nchar(seq.content)
seq.len#输出各组字符数目
a_count<-""
t_count<-""
c_count<-""
g_count<-""
n_count<-""
gc_percent<-""
n_percent<-""#新建变量
for(i in 1:length(seq.content)){
  a_count[i]<-length(gregexpr("[Aa]",seq.content,perl=T)[[i]])
  t_count[i]<-length(gregexpr("[Tt]",seq.content,perl=T)[[i]])
  c_count[i]<-length(gregexpr("[Cc]",seq.content,perl=T)[[i]])
  g_count[i]<-length(gregexpr("[Gg]",seq.content,perl=T)[[i]])
  n_count[i]<-length(gregexpr("[Nn]",seq.content,perl=T)[[i]])
  gc_percent[i]<-(as.numeric(g_count[i])+as.numeric(c_count[i]))/as.numeric(seq.len[i])
  n_percent[i]<-as.numeric(n_count[i])/as.numeric(seq.len[i])
}
result<-data.frame(mychr,seq.len,a_count,t_count,c_count,g_count,n_count,gc_percent,n_percent)
result
回复 支持 反对

使用道具 举报

11

主题

52

帖子

280

积分

中级会员

Rank: 3Rank: 3

积分
280
发表于 2017-1-13 11:36:44 | 显示全部楼层
sxuan 发表于 2017-1-9 21:46
这个upper用后面有啥区别吗

染色体号也被改了,虽然对结果没有太大影响但逻辑是不对的
回复 支持 反对

使用道具 举报

0

主题

6

帖子

291

积分

中级会员

Rank: 3Rank: 3

积分
291
发表于 2017-1-13 13:18:58 | 显示全部楼层
本帖最后由 雨林木风 于 2017-1-13 13:24 编辑

Python-031

解题思路:
1. 逐行读取文件,构建字典以染色体号作为keys,染色体序列作为value存储每一条染色体的数据;
2. 利用len()函数算序列总长,count计算某一碱基在染色体中的数目;
3. 以每条染色体N碱基数目除以总长即为N的比例,G+C的碱基总数处理总长减去N碱基的数目即为G+C的含量。

代码记录:(由于hg19文件太大,个人电脑运行不了,故采用测试文件进行计算)
[Python] 纯文本查看 复制代码
   #! usr/bin/env python
   
   from collections import OrderedDict
   chr = OrderedDict()
   with open('/ifs1/Group2/stu050/python/test.fa') as f:
       for line in f:
           line = line.strip()
           line = line.upper()
           if line.startswith('>'):
              chrname = line
              chr[chrname]=''
          else:
              chr[chrname]+=line
      for seqname,seq in chr.items():
          length = len(seq)
          N =float(seq.count('N'))
          C =float(seq.count('C'))
          G =float(seq.count('G'))
          print seqname,length,"%.2f"%(N/length),"%.2f"%((G+C)/(length-N))


输出结果:
>CHR_1 44 0.09 0.42
>CHR_2 34 0.12 0.43
>CHR_3 33 0.12 0.45
>CHR_4 25 0.12 0.41

小计:由于计算出来的碱基个数和总长都是整型,在算比例时要通过float()函数将其转化为浮点型或者通过乘以1.0的方式将其转换为浮点型,对于该题目后续继续思考尝试更多种途径的解决方案。
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-11-13 03:55 , Processed in 0.046200 second(s), 22 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.