搜索
楼主: Jimmy

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

  [复制链接]

0

主题

6

帖子

53

积分

注册会员

Rank: 2

积分
53
发表于 2017-2-12 16:14:16 | 显示全部楼层
122-R-江  
理解视频上的方法,用的楼主的测试数据
seq<-readLines("hg19_demo.txt")
is.anno<-regexpr("^>",seq,perl = T)#判断数据属于哪一类型
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#start[2:4],但是缺少start[1]
end<-c(end,length(seq))#加上start[1]
distance<-end-start
index<-1:length(start)
index<-rep(index,distance)#以distance作为法则来写,index,2,2,2,1正是其重复次数
seqs<-tapply(seq.content,index,paste,collapse="")#这里用的是collapse,collapse有一组数内部之间关系的意思,而sep好像更多带有单个数后面所添加的东西的含义。
#试了下sep,老是出错,可见在这里使用不对。tapply是对分组数据进行处理。
seq.content<-as.character(seqs)
seq.len<-nchar(seq.content)
A_count<-""
T_count<-""
C_count<-""
G_count<-""
N_count<-""
GC_percent<-""
N_percent<-""#定义,这个应该叫全局变量吧,因为后面for循环要用到它的[i]
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(seq.anno,seq.len,A_count,T_count,C_count,G_count,N_count,GC_percent,N_percent)
result
回复 支持 反对

使用道具 举报

0

主题

3

帖子

99

积分

注册会员

Rank: 2

积分
99
发表于 2017-2-14 15:29:02 | 显示全部楼层
ID:179
Python
第二题确实学到了不少,思路一就是运用字典的方法来进行存储数据,然后提取字典的键和值来进行分析,但会很耗时间。另外一个就是用pysam包。为了解决第一题,电脑安装了64位python3版本的anaconda,使用conda list命令后,发现并没有pysam包。之后使用pip install pysam的命令安装pysam,结果也出现了安装到python2.7的路径下。用vim ~/.bashrc打开文件,在python2.7的路径前加上#后保存退出,再使用source ~/.bashrc进行重置。最后再使用pip install pysam就成功了。最后pysam安装在anaconda3/bin/目录下。pysam确实很好用。
[Python] 纯文本查看 复制代码
import pysam

hg19 = pysam.Fastafile('hg19.fa')
for SeqName in hg19.references:
    seq = hg19.fetch(SeqName)
    seqlength = len(seq)
    N = seq.count('N')
    GC =seq.count('G')+seq.count('g')+seq.count('C')+seq.count('c')
    print(SeqName,seqlength,'%.2f'%(N/seqlength),"%.2f"%(GC/(seqlength-N)))
回复 支持 1 反对 0

使用道具 举报

5

主题

32

帖子

484

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
484
发表于 2017-2-14 22:32:44 | 显示全部楼层
python-167import sys
import time
start=time.clock()
args=sys.argv
sum_atgc={}
f=open('D:/cs.txt')
for line in f:
    if line.startswith('>'):
        chr_id=line[1:]
        sum_atgc[chr_id]={}
        sum_atgc[chr_id]['A']=0
        sum_atgc[chr_id]['T']=0
        sum_atgc[chr_id]['C']=0
        sum_atgc[chr_id]['G']=0
        sum_atgc[chr_id]['N']=0
    else:
        line=line.upper()
        sum_atgc[chr_id]['A']+=line.count('A')
        sum_atgc[chr_id]['T']+=line.count('T')
        sum_atgc[chr_id]['C']+=line.count('C')
        sum_atgc[chr_id]['G']+=line.count('G')
        sum_atgc[chr_id]['N']+=line.count('N')
for chr_id,atgc_count in sum_atgc.items():
    GC=atgc_count['G']+atgc_count['C']
    SUM=sum(atgc_count.values())
    print (chr_id)
    print ('A : %s' % (atgc_count['A']))
    print ('T : %s' % (atgc_count['T']))
    print ('C : %s' % (atgc_count['C']))
    print ('G : %s' % (atgc_count['G']))
    print ('N : %s' % (atgc_count['N']))
    print ('SUM : %s' % (SUM))
    print ('GC : %f' % (GC*1.0/SUM))
    print ('N : %f' % (atgc_count['N']*1.0/SUM))
end=time.clock()   
print('used %s s' % str(end-start))


回复 支持 反对

使用道具 举报

2

主题

9

帖子

235

积分

中级会员

Rank: 3Rank: 3

积分
235
发表于 2017-2-18 17:53:41 | 显示全部楼层
本帖最后由 wangfangce 于 2017-2-18 18:09 编辑

206-perl-策
为了赶上进度,只能直接看jimmy前辈的代码了。。
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl
while (<>){
chomp;
if (/^>/){
$chr=$_;
#print $_;
}else{
$count_A{$chr}+=($_=~tr/Aa//);
$count_C{$chr}+=($_=~tr/Cc//);
$count_G{$chr}+=($_=~tr/Gg//);
$count_T{$chr}+=($_=~tr/Tt//);
$count_N{$chr}+=($_=~tr/Nn//);
}
}
foreach $chr(sort keys %count_N){

$length=$count_A{$chr}+$count_T{$chr}+$count_C{$chr}+$count_G{$chr}+$count_N{$chr};
$GC=($count_C{$chr}+$count_G{$chr})/($count_A{$chr}+$count_T{$chr}+$count_C{$chr}+$count_G{$chr});
print "$chr\t$count_A{$chr}\t$count_T{$chr}\t$count_C{$chr}\t$count_G{$chr}\t$length\t$GC\n";
}

运行结果(部分):
>chr1        65570891        65668756        47024412        47016562        249250621        0.417439252353623
>chr10        38330752        38376915        27308648        27298423        135534747        0.415848760251115
>chr11        38307244        38317436        27236798        27268038        135006516        0.415656502537537
>chr11_gl000202_random        9226        8978        11254        10645        40103        0.546068872652919
>chr12        38604831        38624517        26634995        26617050        133851895        0.40811983820559
>chr13        29336945        29425459        18412698        18414776        115169878        0.385265414817247
>chr14        25992966        26197495        18027132        18071947        107349540        0.408871526570418
>chr15        23620876        23597921        17247582        17228387        102531392        0.42200952016926
>chr16        21724083        21828642        17630040        17701988        903547530.447894259109869
>chr17        21159933        21206981        17727956        17700340        811952100.45540459367614
回复 支持 反对

使用道具 举报

633

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
 楼主| 发表于 2017-2-18 20:33:59 | 显示全部楼层
wangfangce 发表于 2017-2-18 17:53
206-perl-策
为了赶上进度,只能直接看jimmy前辈的代码了。。
[mw_shl_code=perl,true]#!/usr/bin/perl

抓紧时间看基本语法,可以做一些简单的题目,比如:http://www.biotrainee.com/thread-834-1-1.html  ,我很看好你的!
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

1

主题

15

帖子

126

积分

注册会员

Rank: 2

积分
126
发表于 2017-2-21 15:06:17 | 显示全部楼层
本帖最后由 starhqy2 于 2017-2-21 15:08 编辑

199

R:

[Bash shell] 纯文本查看 复制代码
rm(list=ls())

a <- readLines('G:/BaiduYunDownload/第二题/hg19_sub.txt')
pre_seq <- strsplit(paste0(sub('(^>.*)', ':\\1,', a), collapse = ''), ':')[[1]]#先添加几个分割点, 然后合并所有字符,最后根据序列再分割
s <- strsplit(pre_seq, ',') # s是包含每条序列的列表
output <- data.frame(chr = '', A = '', G = '', C = '', T = '', N = '', GC = '', N_percent = '') #创建数据框用于存储每条序列的ATGC信息
lapply(s, function(line){
  chr <- line[1]
  A <- length(gregexpr('[Aa]', line[2])[[1]])
  T <- length(gregexpr('[Tt]', line[2])[[1]])
  C <- length(gregexpr('[Cc]', line[2])[[1]])
  G <- length(gregexpr('[Gg]', line[2])[[1]])
  N <- length(gregexpr('[Nn]', line[2])[[1]])
  GC <- (as.numeric(G) + as.numeric(C))/as.numeric(length(line[2]))
  N_percent <- as.numeric(N)/as.numeric(length(line[2]))
  output <<- rbind(output, data.frame(chr = line[1], A=as.character(A), G = as.character(G), C = as.character(C),  T= as.character(T), N= as.character(N) ,  GC = as.character(GC) , N_percent = as.character(N_percent))) #使用全局定义变量, 合并每次分析的结果
})output


回复 支持 反对

使用道具 举报

0

主题

3

帖子

131

积分

注册会员

Rank: 2

积分
131
发表于 2017-2-25 00:17:03 | 显示全部楼层
136-R+python-Yangk
拖得有点久,逐渐赶回来
seq<-readLines('G:/biotrainee/2nd/hg19.txt')
is.title<-regexpr('^>',seq,perl=T)
seq.name<-which(is.title==1)
seq.content<-which(is.title==-1)
start<-seq.name
end<-start[2:length(start)]-1
end<-c(end,length(seq))
distance<-end-start
index<-rep(1:length(start),distance)
seq2<-tapply(seq[seq.content],index,paste,collapse='')
seq2<-as.character(seq2)
seq_length<-nchar(seq2)
A_count<-""
T_count<-""
C_count<-""
G_count<-""
GC_percent<-""
for(i in 1:length(start))
{A_count[i]=length(gregexpr('[Aa]',seq2,perl = T)[[i]])
T_count[i]=length(gregexpr('[Tt]',seq2,perl = T)[[i]])
C_count[i]=length(gregexpr('[Cc]',seq2,perl = T)[[i]])
G_count[i]=length(gregexpr('[Gg]',seq2,perl = T)[[i]])
GC_percent[i]=(as.numeric(G_count[i])+as.numeric(C_count[i]))/as.numeric(seq_length[i])
}
results<-data.frame(seq[seq.name],seq_length,A_count,T_count,G_count,C_count,GC_percent)
results


运行结果
seq.seq.name. seq_length A_count T_count G_count C_count        GC_percent
1        >chr_1         44      13      10      10       7 0.386363636363636
2        >chr_2         34      11       6       8       5 0.382352941176471
3        >chr_3         33      10       6      10       3 0.393939393939394
4        >chr_4         25       9       4       7       2              0.36
回复 支持 反对

使用道具 举报

0

主题

3

帖子

208

积分

中级会员

Rank: 3Rank: 3

积分
208
发表于 2017-2-26 16:28:48 | 显示全部楼层
我用了biopython模块,处理一些文件也挺快的
[Python] 纯文本查看 复制代码
from __future__ import division
import sys, getopt
import operator
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
import time

start = time.clock()

out_file = open('hg19_result.txt','w+')
for record in SeqIO.parse("hg19.fasta","fasta"):
    seq = record.seq
    Len = len(seq)
    N = seq.count("N")
    GC = seq.count("G") + seq.count("g") + seq.count("C") + seq.count("c")
    print Len, N, GC
    result = "%s%s%s%s%.2f%s%.2f%s" % (record.id,'\t',Len,'\t',N/Len,'\t',GC/(Len-N),'\n')
    out_file.write(result)
end = time.clock()
out_file.write("The function run time is : %.03f seconds" %(end-start))

   
out_file.close

回复 支持 反对

使用道具 举报

0

主题

3

帖子

208

积分

中级会员

Rank: 3Rank: 3

积分
208
发表于 2017-2-26 16:36:20 | 显示全部楼层
142号学员,昨天的也忘了加了
回复 支持 反对

使用道具 举报

5

主题

32

帖子

484

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
484
发表于 2017-2-26 21:52:20 | 显示全部楼层
#! /usr/bin/perl
use strict;
open INPUT,"chr1.fa";
open OUT,">result";
my (@gene,$chr,%A_num,%T_num,%G_num,%C_num,%N_num,$len,$gc,$n);
chomp(@gene=<INPUT>);
foreach (@gene){
if(/>/){
$chr=$_;
}
else{
$A_num{$chr} +=($_=~tr/Aa//);
$T_num{$chr} +=($_=~tr/Tt//);
$C_num{$chr} +=($_=~tr/Cc//);
$G_num{$chr} +=($_=~tr/Gg//);
$N_num{$chr} +=($_=~tr/Nn//);
}
}
foreach(sort keys %A_num){
$len =$A_num{$_}+$T_num{$_}+$G_num{$_}+$C_num{$_}+$N_num{$_};
$gc=sprintf "%.2f",(($G_num{$_}+$C_num{$_})/$len);
$n=sprintf "%.2f",($N_num{$_}/$len);
print OUT "$_\t$A_num{$_}\t$T_num{$_}\t$G_num{$_}\t$C_num{$_}\t$N_num{$_}\t$len\t$gc\t$n\n";
}


#####
自己电脑不太好,只计算了第一个染色体的数量
>chr1        65570891        65668756        47016562        47024412        23970000        249250621        0.38        0.10
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-12-8 14:40 , Processed in 0.029666 second(s), 22 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.