搜索
楼主: Jimmy

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

  [复制链接]

0

主题

3

帖子

48

积分

新手上路

Rank: 1

积分
48
发表于 2017-1-17 19:53:38 | 显示全部楼层
本帖最后由 perlyi 于 2017-1-17 19:54 编辑

学员编号:097

python

from __future__ import division
from collections import OrderedDict
dict = OrderedDict()
name = ""
with open('C:\\1\\1.txt') as f:
    for line in f:
        context = ""
        line=line.strip()
        if line.startswith(">"):
            name = line
            dict[name]=""
        else:context = line
        dict[name]=dict[name]+context
for name in dict:
    print name,dict[name]
    print "length=",len(dict[name])
    print "N is",round(dict[name].count("N")/len(dict[name]),2)
    print "G is",round((dict[name].count("G")+dict[name].count("g"))/len(dict[name]),2)
    print "C is",round((dict[name].count("C") + dict[name].count("c")) / len(dict[name]),2)

回复 支持 反对

使用道具 举报

0

主题

3

帖子

46

积分

新手上路

Rank: 1

积分
46
发表于 2017-1-18 14:45:49 | 显示全部楼层
082-Python/R-wiwi

[mw_shl_code=python,true]#第一种方法
with open ("myfastq.fasta","r") as hg19:
    for line in hg19:
        line=line.strip()
        if line.startwith(">"):
            tmp_chr=line
            chr_dict[tmp_chr]=""
        else:
            chr_dict[tmp_chr]+=line

for seqName, seq in chr_dict.items():
    seqLen=len(seq)
    N=seq.count("N")
    GC=seq.count("G")+seq.count("g")+seq.count("C")+seq.count("c")
    print(seqName,seqLen,"%.2f"%(N/seqLen),"%.2f"%(GC/seqLen-N))
[/mw_shl_code]
回复 支持 反对

使用道具 举报

13

主题

26

帖子

344

积分

中级会员

Rank: 3Rank: 3

积分
344
发表于 2017-1-19 10:10:06 | 显示全部楼层
045-python-R
前面学的又忘了,对循环还是理解的不透彻啊,先用测试数据跑一遍,回头再把基础知识看一遍[mw_shl_code=python,true]import os
from collections import OrderedDict      #collections模块带有OrderedDict有序字典
os.chdir("E:/python")                               #改变到当前工作目录

chr_dict = OrderedDict()                         #定义chr_dict做为有序字典
temp_dict = OrderedDict()                     #定义temp_dict做为有序字典
temp_chr = ""

with open("test.fasta","r") as test:           #with打开文件
    for line in test:
        line = line.strip()                               #strip()函数去除每行首位的空字符
        if line.startswith(">"):                       #startswith()函数判断行开头是否是‘>’
            temp_chr = line                           #如果是‘>’,把这行赋值给temp_chr
            chr_dict[temp_chr] = ""               
        else:
            chr_dict[temp_chr] += line
#print(chr_dict)                                         #绕糊涂了,打印出来看看
#OrderedDict([('>chr_1', 'ATCGTCGaaAATGAANccNNttGTAAGGTCTNAAccAAttGggG'), ('>chr_2', 'ATCGAATGATCGANNNGccTAAGGTCTNAAAAGG'), ('>chr_3', 'ATCGTCGANNNGTAATggGAAGGTCTNAAAAGG'), ('>chr_4', 'ATCGTCaaaGANNAATGANGgggTA')])
for seqName, seq in chr_dict.items():
    seqLen = len(seq)
    seq = seq.upper()                                #将seq序列全部转换成大写
    N = seq.count("N")
    GC = seq.count("G") + seq.count("C")
    print(seqName, seqLen, "%.2f"%(N/seqLen), "%.2f"%(GC/(seqLen-N)))  #定义输出结果保留两位小数[/mw_shl_code]
结果:>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
回复 支持 反对

使用道具 举报

1

主题

15

帖子

180

积分

注册会员

Rank: 2

积分
180
发表于 2017-1-20 19:43:23 | 显示全部楼层
AdaWon 发表于 2017-1-16 21:32
请教一下,这两行代码什么意思?
$/ = ">";;
my($chr,$seq) = split;

我设置了起读为> 就是每次读取不是按换行符(默认是换行符),而是从这个>读到下一次>前。
举个例子,假设序列是这样的:

>23
aaagggcccgggg
>3455
aaggggggg
那么设置$/ = ">"第一次读到的内容就是:
>23
aaagggcccaaggggg

my($chr,$seq) = split;的意思是split默认分割,分割成的两部分分别是$chr,$seq

回复 支持 反对

使用道具 举报

0

主题

19

帖子

241

积分

中级会员

Rank: 3Rank: 3

积分
241
发表于 2017-1-21 22:56:08 | 显示全部楼层
本帖最后由 朝曦曦 于 2017-1-22 17:54 编辑

学员205
perl
首先跟着视频完成测试文件。命名为test_2.txt,print看是否正确 perl -lane '{print}' test_2.fa写出抓取ATCG数量并计算GC百分含量的脚本(模仿群主写的,以前有的不理解,经过视频讲解,几乎都理解了)命名为lesson2_1.pl
[mw_shl_code=applescript,true]while (<>){                     #reading file
chomp;                          #delete white space

if(/^>/){
$chr=$_                        #if meet >beginning line, save it to $chr
}
else{
$A_count{$chr}+=($_=~tr/Aa//); #counting how many Aa and save the result to $A_count
$T_count{$chr}+=($_=~tr/Tt//);
$C_count{$chr}+=($_=~tr/Cc//);
$G_count{$chr}+=($_=~tr/Gg//);
$N_count{$chr}+=($_=~tr/Nn//);
}
}
foreach (sort keys  %N_count){
$length = $A_count{$_}+$T_count{$_}+$C_count{$_}+$G_count{$_}+$N_count{$_};
$gc = ($G_count{$_}+$C_count{$_})/($A_count{$_}+$T_count{$_}+$C_count{$_}+$G_count{$_});
print "$_\t$A_count{$_}\t$T_count{$_}\t$C_count{$_}\t$G_count{$_}\t$N_count{$_}\t$length\t$gc\n"

}[/mw_shl_code]
运行perl lesson2_1.pl test_2.txt 看是否成功,并计算正确。确认无误后,再运行基因组计算。
我没有能够运行成功网页中给的下载hg19的代码,在网上搜了一个,也能下载就没再改。
[mw_shl_code=applescript,true]#download hg19 as the ftp://hgdownload.cse.ucsc.edu/go ... omosomes/README.txt website instruction
wget --timestamping
        'ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/*'[/mw_shl_code]

最后得到运算结果
perl lesson2_1.pl hg19/chr1.fa
>chr1        65570891        65668756        47024412        47016562        2397000249250621        0.417439252353623
chr1为例。


回复 支持 反对

使用道具 举报

0

主题

4

帖子

51

积分

注册会员

Rank: 2

积分
51
发表于 2017-1-23 17:58:01 | 显示全部楼层
本帖最后由 zhongzheng 于 2017-1-23 18:13 编辑

212-perl-zhongzheng
视频中的代码简洁明了,下面的代码稍微繁琐了些:

[mw_shl_code=perl,true]#!/usr/bin/perl
use warnings;
use strict;

my ($row,$chr,%A_count,%T_count,%G_count,%C_count,%N_count,$length,$gc);
open DATA, "h19.fa";
open OUTPUT, ">atcg.txt";

while ($row=<DATA>) {
chomp ($row);
    last unless $row =~ /\S/;
    if ($row =~ /^>/){
        $chr = $row;
    }
    $A_count{$chr} += $row=~tr/Aa//;
    $T_count{$chr} += $row=~tr/Tt//;
    $G_count{$chr} += $row=~tr/Gg//;
    $C_count{$chr} += $row=~tr/Cc//;
    $N_count{$chr} += $row=~tr/Nn//;        
}

print OUTPUT "chrome\tA\tT\tG\tC\tN\tLength\tGC\%\n";
  
foreach $chr(sort keys %A_count){
    $length=$A_count{$chr}+$T_count{$chr}+$G_count{$chr}+$C_count{$chr}+$N_count{$chr};
        $gc=($G_count{$chr}+$C_count{$chr})*100/($A_count{$chr}+$T_count{$chr}+$G_count{$chr}+$C_count{$chr});
    print OUTPUT "$chr\t$A_count{$chr}\t$T_count{$chr}\t$G_count{$chr}\t$C_count{$chr}\t$N_count{$chr}\t$length\t$gc\n";
}[/mw_shl_code]
运行结果:
hg19.png
回复 支持 反对

使用道具 举报

0

主题

9

帖子

79

积分

注册会员

Rank: 2

积分
79
发表于 2017-1-24 18:43:57 | 显示全部楼层
010-python
模仿之前有一个练习的作业,计算那个基因的GC含量多写的
[mw_shl_code=python,true]import re
import os
from collections import OrderedDict
from operator import itemgetter
os.chdir('C:/Users/Administrator/Desktop/')
seq=OrderedDict()
with open ('test fastq data.txt') as f:
    for line in f:
        line=line.rstrip()
        if line.startswith('>'):
            seqName=line[1:]
            seq[seqName]=''
        else:
            seq[seqName]+=line.upper()
for key,value in seq.items():
    seqlength=len(value)
    Nnumber=value.count('N')
    GCnumber=value.count('G')+value.count('C')
    GCcontent=float(GCnumber/seqlength)*100
    print(seqname,'\t',seqlength,'\t',Nnumber,'\t',GCcontent)[/mw_shl_code]
chr_1          44          4          38.63636363636363
chr_2          34          4          38.23529411764706
chr_3          33          4          39.39393939393939
chr_4          25          3          36.0
后来看了视频,感觉用if 。。。else,更加的简洁明了
回复 支持 反对

使用道具 举报

0

主题

6

帖子

131

积分

注册会员

Rank: 2

积分
131
发表于 2017-1-30 01:01:51 | 显示全部楼层
074-python

下载了示例demo数据,用视频中老师提到的方法进行了学习
1.普通思路[mw_shl_code=python,true]#!/usr/bin/python
#-*-coding:utf-8-*-


from collections import OrderedDict
chr_dict = OrderedDict()
tem_chr = ''

with open('hg19_demo.txt','r') as hg19:
    for line in hg19:
        line=line.strip()
        if line.startswith('>'):
            tem_chr = line   #如果line以>开头,把该line写进字典的键
            chr_dict[tem_chr]= '' #给字典的键一个初始值
        else:
            chr_dict[tem_chr] += line
for seqName, seq in chr_dict.items(): #dict.items()以列表返回可遍历的(键, 值) 元组数组
    seqLen = len(seq)
    N = float(seq.count('N'))       #下面涉及保留两位小数,这里用float函数限定数据为浮点格式
    G = float(seq.count('G') + seq.count('g'))
    C = float(seq.count('C') + seq.count('c'))
    print(seqName,seqLen,"%.2f"%(N/seqLen),"%.2f"%((G+C)/(seqLen-N)))  #%.2f 浮点数格式,保留两位小数。中间蓝色的%表示格式化操作[/mw_shl_code]
2.用pysam模块进行分析
[mw_shl_code=applescript,true]#!/usr/bin/python
#-*-coding:utf-8-*-

import pysam


hg19 = pysam.FastaFile("hg19_demo.fasta")
dir(hg19) #list methods
list(zip(hg19.references,hg19.lengths))
for seqName in hg19.references:
    seq = hg19.fetch(seqName)
    seqLen = len(seq)
    N = float(seq.count('N'))      
    G = float(seq.count('G') + seq.count('g'))
    C = float(seq.count('C') + seq.count('c'))
    print(seqName,seqLen,"%.2f"%(N/seqLen),"%.2f"%((G+C)/(seqLen-N))) [/mw_shl_code]
结果:
('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')



回复 支持 反对

使用道具 举报

0

主题

6

帖子

131

积分

注册会员

Rank: 2

积分
131
发表于 2017-1-30 01:16:55 | 显示全部楼层
熊海清 发表于 2017-1-30 01:01
074-python

下载了示例demo数据,用视频中老师提到的方法进行了学习

类似方法得到人基因组结果如下,运行非常快
('chrM', 16571, '0.00', '0.44')
('chr1', 249250621, '0.10', '0.42')
('chr2', 243199373, '0.02', '0.40')
('chr3', 198022430, '0.02', '0.40')
('chr4', 191154276, '0.02', '0.38')
('chr5', 180915260, '0.02', '0.40')
('chr6', 171115067, '0.02', '0.40')
('chr7', 159138663, '0.02', '0.41')
('chr8', 146364022, '0.02', '0.40')
('chr9', 141213431, '0.15', '0.41')
('chr10', 135534747, '0.03', '0.42')
('chr11', 135006516, '0.03', '0.42')
('chr12', 133851895, '0.03', '0.41')
('chr13', 115169878, '0.17', '0.39')
('chr14', 107349540, '0.18', '0.41')
('chr15', 102531392, '0.20', '0.42')
('chr16', 90354753, '0.13', '0.45')
('chr17', 81195210, '0.04', '0.46')
('chr18', 78077248, '0.04', '0.40')
('chr19', 59128983, '0.06', '0.48')
('chr20', 63025520, '0.06', '0.44')
('chr21', 48129895, '0.27', '0.41')
('chr22', 51304566, '0.32', '0.48')
('chrX', 155270560, '0.03', '0.39')
('chrY', 59373566, '0.57', '0.40')
回复 支持 1 反对 0

使用道具 举报

2

主题

21

帖子

245

积分

中级会员

Rank: 3Rank: 3

积分
245
发表于 2017-2-2 19:51:36 | 显示全部楼层
disheng 发表于 2017-1-20 19:43
我设置了起读为> 就是每次读取不是按换行符(默认是换行符),而是从这个>读到下一次>前。
举个例子,假 ...

谢谢你!
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2020-2-29 08:10 , Processed in 0.025180 second(s), 22 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.