搜索
楼主: Jimmy

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

  [复制链接]

0

主题

22

帖子

183

积分

注册会员

Rank: 2

积分
183
发表于 2017-1-13 16:01:05 | 显示全部楼层
okwell 发表于 2017-1-11 21:47
052-perl+R
才开始接触生信编程,懂得的还很少,只好根据论坛上的答案来敲了,有点抄作业的感觉
思路:计 ...

加油!代码写不出没关系,可以先写一段自己的思路总结,然后再给教学视频的代码写注解,最后自己再尝试能不能写出来。
回复 支持 反对

使用道具 举报

0

主题

2

帖子

22

积分

新手上路

Rank: 1

积分
22
发表于 2017-1-14 00:33:16 | 显示全部楼层
刚填了第一题的坑(虽然没填完)
顺便爬来第二题再来挖新坑(莫名觉得自己有好多坑没填???)
照惯例陈述原因+立个flag
期末考到17号然后会留校做一周实验再回家,争取23号前,填完这个坑!
回复 支持 反对

使用道具 举报

0

主题

2

帖子

176

积分

注册会员

Rank: 2

积分
176
发表于 2017-1-14 03:11:53 | 显示全部楼层
本帖最后由 juzheng87 于 2017-1-14 03:15 编辑

148 R-shell-perl-python
正在学习中,这一周刚把第一题视频中R,perl,python的内容整明白,感觉perl处理这个问题更顺畅。正在进行模仿。争取25号之前完成第一题和第二题的内容,赶上进度。

回复 支持 反对

使用道具 举报

1

主题

15

帖子

126

积分

注册会员

Rank: 2

积分
126
发表于 2017-1-14 10:36:57 | 显示全部楼层
本帖最后由 starhqy2 于 2017-1-14 10:54 编辑

不用字典, 用字典太占内存了。by: Padre
[Python] 纯文本查看 复制代码
#!/usr/bin/env python
#-*- coding:utf-8 -*-

############## Usage #########
#
#      count_atgc.py hg19.fa
#
###################################

import sys
chr_name = ''
a = 0
t = 0
c = 0
g = 0
gc = 0
length = 0
flag = False

with open('ATCG_count.txt', 'w') as f_save:
    print('{0}\t{1}\t{2}\t{3}\t{4}\t{5}'.format('chr_name', 'A', 'T', 'C', 'G', 'GC_percent'), file = f_save)
    with open(sys.argv[1], 'r') as f_genome:
        for line in f_genome:
            line = line.strip()
            if line.startswith('>'):
                if flag:
                    print('{0}\t{1}\t{2}\t{3}\t{4}\t{5}'.format(chr_name, a, t, c, g, round(gc/length,4)), file = f_save)
                    chr_name = ''
                    a = 0
                    t = 0
                    c = 0
                    g = 0
                    gc = 0
                    length = 0
                chr_name = line[1:]
            else:
                flag = True
                line = line.upper()
                a += line.count('A')
                t += line.count('T')
                c += line.count('C')
                g += line.count('G')
                gc += g + c
                length += len(line)
        if flag:
            print('{0}\t{1}\t{2}\t{3}\t{4}\t{5}'.format(chr_name, a, t, c, g, round(gc/length,4)), file = f_save)

回复 支持 反对

使用道具 举报

4

主题

14

帖子

105

积分

注册会员

Rank: 2

积分
105
发表于 2017-1-14 10:52:52 | 显示全部楼层
from collections import OrderedDict
chr_dict={}
temp_chr = " "
with open('no_2.txt','r') as no_2:
    for line in no_2:
        line = line.strip()
        if line.startswith(">"):
            temp_chr = line
            chr_dict[temp_chr] = ""
        else:
            chr_dict[temp_chr] += line
for seqname, seq in chr_dict.items():
    print(seqname,seq)
    seqLen = len(seq)
    N  =  seq.count('N')
    GC = seq.count('G') +seq.count('C')+seq.count('g')+seq.count('c')
    print(seqLen,N.GC/(seqLen-N))

    看了视频之后,发现思路都差不多,打开文件,读取,选择符合条件的行,然后对其进行操作,行名称与行内容一一对应,统计行内碱基个数,按行号输出统计信息等。
现在学会了处理的思路,但是如何用python实现这些要求还有些力不能及,比如行号与行对应添加到字典中,自己独立做还做不出。
    另外,用测试数据可以进行处理,对下载来的基因组数据同样的代码却出现错误,chr_dict[temp] += line报错Keyerror = ‘’,这一点还不能理解。
    除此之外,还学会了几个linux命令:
wget + url 下载
tar:压缩及解压缩。-vczf各种参数
rm:删除文件   -rf 强制删除
*.fa:通配符*
cat:读取文件内容
回复 支持 反对

使用道具 举报

0

主题

6

帖子

46

积分

新手上路

Rank: 1

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

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

终于下载了hg19的数据,但是里面不同的染色体分成了不同的文件,所以,在读取文件的时候需要读取文件夹中的文件,要import os,

在学习了dongye老师的教学视频时候,修改了代码,没有像之前的代码一样将碱基信息放入字典中,而是识别ATCG直接计数,这样计算结果为:

chr1 249250621
N: 0.1
GC: 0.38
AT: 0.49

chr10 135534747
N: 0.03
GC: 0.4
AT: 0.53

chr11 135006516
N: 0.03
GC: 0.4
AT: 0.52

chr12 133851895
N: 0.03
GC: 0.4
AT: 0.54

chr13 115169878
N: 0.17
GC: 0.32
AT: 0.47

chr14 107349540
N: 0.18
GC: 0.34
AT: 0.45

chr15 102531392
N: 0.2
GC: 0.34
AT: 0.43

chr16 90354753
N: 0.13
GC: 0.39
AT: 0.46

chr17 81195210
N: 0.04
GC: 0.44
AT: 0.5

chr18 78077248
N: 0.04
GC: 0.38
AT: 0.53

chr19 59128983
N: 0.06
GC: 0.46
AT: 0.47

chr2 243199373
N: 0.02
GC: 0.39
AT: 0.54

chr20 63025520
N: 0.06
GC: 0.42
AT: 0.5

chr21 48129895
N: 0.27
GC: 0.3
AT: 0.4

chr22 51304566
N: 0.32
GC: 0.33
AT: 0.34

chr3 198022430
N: 0.02
GC: 0.39
AT: 0.55

chr4 191154276
N: 0.02
GC: 0.38
AT: 0.56

chr5 180915260
N: 0.02
GC: 0.39
AT: 0.55

chr6 171115067
N: 0.02
GC: 0.39
AT: 0.55

chr7 159138663
N: 0.02
GC: 0.4
AT: 0.54

chr8 146364022
N: 0.02
GC: 0.39
AT: 0.54

chr9 141213431
N: 0.15
GC: 0.35
AT: 0.46

chrM 16571
N: 0.0
GC: 0.44
AT: 0.56

chrX 155270560
N: 0.03
GC: 0.38
AT: 0.53

chrY 59373566
N: 0.57
GC: 0.17
AT: 0.23

计时:
real        6m19.964s
user        6m6.884s
sys        0m3.009s

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

filepath = '/Users/sunshine/practice/practice2-hg19/chromFa/'
fileDict = []
filenames = os.listdir(filepath)

for f in filenames:
    if f.startswith('.'): #remove hidden files
        pass
    else:
        fileDict.append(f)

chr_dict = collections.OrderedDict()
bases = ['A','T','C','G','N','a','t','c','g','n']

for f in fileDict:
    f = open(filepath+f,'r')
    for line in f:
        line = line.strip('\n')
        if line.startswith('>'):
            temp = line[1:]
            chr_dict[temp] = {}
            for base in bases:
                chr_dict[temp][base] = 0
        else:
            for base in bases:
                chr_dict[temp][base] += line.count(base)
                
for chr_name, atcgn in chr_dict.items():
    length = sum(atcgn.values())
    N = round((atcgn['N']+atcgn['n'])/length,2)
    GC = round((atcgn['G']+atcgn['C']+atcgn['g']+atcgn['c'])/length,2)
    AT= round((atcgn['A']+atcgn['T']+atcgn['a']+atcgn['c'])/length,2)
    print(chr_name, length)
    print('N: '+ str(N))
    print('GC: '+ str(GC))
    print('AT: '+ str(AT))
    print('')
    
f.close()
回复 支持 3 反对 0

使用道具 举报

6

主题

36

帖子

465

积分

中级会员

Rank: 3Rank: 3

积分
465
发表于 2017-1-14 16:01:06 | 显示全部楼层
本帖最后由 017中大普外 于 2017-1-14 16:02 编辑

感觉还是对字典不熟,一些灵活的用法还需要掌握
[Python] 纯文本查看 复制代码
import os
[mw_shl_code=python,true]import os
import re

os.chdir("F:\生信菜鸟团作业")
f1=open("dierci.txt","r")

for i in f1.readlines():
    i=i.strip("\n")
    chroma=re.split(">",i)
    print(chroma)
    print(i)
    if i.startswith(">"):
        print()
        
with open ('F:\生信菜鸟团作业\dierci.txt') as f:
    for line in f:
        if line.startswith('>'):print(line);continue
        if len(line.strip()) == 0:continue
        length=len(line.strip('\n'))
        print("N:%f\nC+G:%f "%(float(line.count('N'))/length,float(line.count('C')+line.count('G'))/length))
回复 支持 反对

使用道具 举报

0

主题

4

帖子

26

积分

新手上路

Rank: 1

积分
26
发表于 2017-1-14 17:28:17 | 显示全部楼层
036秋天,这周作业交晚了,因为年终工作总结。
while (<>){
chomp;#读取文件,去掉文件内的换行符。

if(/^>/){#调用if函数如果开头是>就进行下面的分析,因为是fasta文件,所以每条序列以>开头,就是导入每条序列。
$chr=$_#将默认标量赋值给chr标量,也就是把染色体的名字赋值给chr
}
else{ #将不是以#开头的数据赋值给以下标量。
$A_count{$chr}+=($_=~tr/Aa//);#对A or a进行统计
$T_count{$chr}+=($_=~tr/Tt//);#对T or t进行统计
$C_count{$chr}+=($_=~tr/Cc//); #对C or c进行统计
$G_count{$chr}+=($_=~tr/Gg//);#对G or g进行统计
$N_count{$chr}+=($_=~tr/Nn//); #对N的未测通序列进行统计
}
}
foreach (sort keys  %N_count){调用foreach函数来对排序后的哈希进行循环
$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#得到GC含量,利用哈希进行统计。
"$_\t$A_count{$_}\t$T_count{$_}\t$C_count{$_}\t$G_count{$_}\t$N_count{$_}\t$length\t$gc\n"
#打印出每一个碱基的数量,统计N值,以分隔符分割。
}
回复 支持 反对

使用道具 举报

0

主题

2

帖子

30

积分

新手上路

Rank: 1

积分
30
发表于 2017-1-14 18:22:44 | 显示全部楼层
本帖最后由 068-R 于 2017-1-14 18:24 编辑

068-R
本题难点:
怎么把数据排列转成能代码能计算的行列

继续学习视频中R语言,并注释
首先在txt文本中输入
>chr_1
ATCGTCGaaAATGAANccNNttGTA
AGGTCTNAAccAAttGggG
>chr_2
ATCGAATGATCGANNNGccTA
AGGTCTNAAAAGG
>chr_3
ATCGTCGANNNGTAATggGA
AGGTCTNAAAAGG
>chr_4
ATCGTCaaaGANNAATGANGgggTA

测试数据,保存为hg_19_demo
[AppleScript] 纯文本查看 复制代码
seq <- readLines("hg_19_demo.txt")  #读取测试数据
is.anno <- regexpr("^>", seq, perl =t)   #用正则表达式匹配有“^>”的染色体号行,染色体号行为1,否者为-1
seq.anno <- seq[which(is.anno == 1)]  #染色体内容
seq.content <- seq[which(is.anno==-1)] #碱基序列内容
start <- which(is.anno == 1) #染色体行号
end <- start[2:length(start)] - 1 #第二条染色体行号到最后条染色体行号-1 即为个染色体中最后序列的行号,但不包含最后一条序列行号
end <- c(end, length(seq))   #末尾加一行,得到每条染色体碱基序列结束行
distance <- end - start  #各染色体中碱基序列所占行数
index <- 1:length(start) #染色体数
index <- rep(index,distance) #重复染色体行数,有多少行碱基对应的染色体数重复多少次
seqs <- tapply(seq.content, index, paste,collapse="") #以index为索引,做分组运算,输出染色体号为head的数组
seq.content <- as.character(seqs) #转变成字符串向量
seq.len <- nchar(seq.content) #计算每条染色体的碱基总数
A_count <- ""
T_count <- ""
C_count <- ""
G_count <- ""
N_count <- "" 
GC_percent <- ""
N_percent <- "" #建立空向量
CA_count <- gregexpr("[Aa]",seq.content,perl = t)
CT_count <- gregexpr("[Tt]",seq.content,perl = t)
CC_count <- gregexpr("[Cc]",seq.content,perl = t)
CG_count <- gregexpr("[Gg]",seq.content,perl = t)
CN_count <- gregexpr("[Nn]",seq.content,perl = t) #匹配向量seq.content中各元素碱基位置,输出为list
for(i in 1: length(seq.content))
{
  A_count[i] <- length(CA_count[[i]])
  T_count[i] <- length(CT_count[[i]])
  C_count[i] <- length(CC_count[[i]])
  G_count[i] <- length(CG_count[[i]])
  N_count[i] <- length(CN_count[[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(seq.anno,seq.len,A_count,T_count,C_count,G_count,N_count,GC_percent,N_percent)

#计算各染色体中碱基的个数及百分比,生成数据框


最后输出
  seq.anno seq.len A_count T_count C_count G_count N_count        GC_percent          N_percent
1   >chr_1      44      13      10       7               4 0.386363636363636 0.0909090909090909
2   >chr_2      34      11       6       5               4 0.382352941176471  0.117647058823529
3   >chr_3      33      10       6       3               4 0.393939393939394  0.121212121212121
4   >chr_4      25       9       4       2               3              0.36               0.12


通过本题训练,掌握了2个正则表达式的用法,学习并回顾apply家族函数循环运算


回复 支持 反对

使用道具 举报

1

主题

5

帖子

87

积分

注册会员

Rank: 2

积分
87
发表于 2017-1-14 19:23:22 | 显示全部楼层
191-perl+R-芃


[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;
my %hash;

open(IN,"$ARGV[0]");

$/=">";

while(<IN>){
    chomp;
    my @a=split("\n",$_);
    foreach my $a (@a){
        if($a=~/chr/){
            print "$a\t";
        }else{
            my $up=uc $a;
            my @str=split("",$up);
            foreach my $b (@str){
                $hash{$b}++;
            }
        }
    }
    foreach my $temp (keys %hash){
        print "$temp:$hash{$temp}\t";
    }
    print "\n";
%hash=();
}

close IN;




hg19文件实在太大了,答案就先不附后面了
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-9-21 01:13 , Processed in 0.041310 second(s), 21 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.