搜索
楼主: Jimmy

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

  [复制链接]

0

主题

2

帖子

85

积分

注册会员

Rank: 2

积分
85
发表于 2017-3-21 18:35:17 | 显示全部楼层
R语言,重新发一下
[mw_shl_code=python,true]
##速度挺快的,可以试试
rm(list = ls())
# 开始计时
time1 <- Sys.time()
#读入数据
data <- read.csv("chrY.fa",header = T,stringsAsFactors = F)
#数据大写
data[,1] <- toupper(data[,1])
#建立函数
get_ATCGN <- function(data_input){
  x <- unlist(strsplit(data_input,split = character(0)))
  A_sum <- sum(x == "A")
  T_sum <- sum(x == "T")
  C_sum <- sum(x == "C")
  G_sum <- sum(x == "G")
  N_sum <- sum(x == "N")
  ALL_sum <- sum(!x == "")
  c(A_sum,T_sum,C_sum,G_sum,N_sum,ALL_sum)
}
#统计结果
x <- get_ATCGN(data[,1])
names(x) <- c("A","T","C","G","N","Total")
x
#计算时间
difftime(Sys.time(),time1, units = "secs")[/mw_shl_code]
回复 支持 反对

使用道具 举报

0

主题

1

帖子

37

积分

新手上路

Rank: 1

积分
37
发表于 2017-6-2 22:54:34 | 显示全部楼层
刚刚开始接触生物信息,跟着老师录制的86分钟的python视频一路学下来收获良多,但是关于第三个部分提高效率的处理有一些疑问,想请教一下。[mw_shl_code=applescript,true]#-*- coding:utf-8 -*-
import re
sum_atgc = {}
bases = ('A','T','C','G')
with open('/Users/user/Desktop/sequence.fasta','r') as Fin:
    tmp = Fin.readline()
    chr_id = re.split(r'\s', tmp)[0][1:]
    sum_atgc[chr_id] = {}#chr_id存入第一条染色体的名称
    for base in bases:
        sum_atgc[chr_id][base] = 0
    print tmp
    print chr_id
    while 1:
        buffer = Fin.read(1024*1024)
        print buffer
        if not buffer:
            break
        if '>' in buffer:#以下处理文件中的第二种情况
            (buffer,tmp) = buffer.split('>')
            print buffer
            print tmp
            buffer = buffer.upper()
            for base in bases:
                sum_atgc[chr_id][base] += buffer.count(base)
                print sum_atgc[chr_id][base]#这里输出的是第一条染色体的碱基数目,与第二条染色体无关
            (tmp,buffer) = tmp.split('\n',1)
            chr_id = re.split(r'\s', tmp)[0]
            print tmp
            print buffer
            print chr_id
            sum_atgc[chr_id] = {}
            for base in bases:
                sum_atgc[chr_id][base] = 0
            buffer = buffer.upper()
            for base in bases:
                sum_atgc[chr_id][base] += buffer.count(base)
                print sum_atgc[chr_id][base]
        else:
            buffer = buffer.upper()
            for base in bases:
                sum_atgc[chr_id][bases] += buffer.count(base)
            print sum_atgc[chr_id][bases][/mw_shl_code]
因为不太理解老师写的程序,所以我用了一个小一些的fasta文件测试了一下,前后有一些代码没有加上,主要的目的是为了弄懂buffer的相关操作,因此基本上print出了每一步输入的结果。但是还是有几个问题不太清楚:

1)if '>' in buffer:
    (buffer,tmp) = buffer.split('>')我一开始用的序列中有两条染色体的序列信息,但是当我改为包含三条序列信息的文件时,因为.split('>')的操作应该是将我的序列信息分成了三个字符串,但是我前面输入的元祖中只有两个元素,因此我的python会提示我ValueError: too many values to unpack,但是看老师的处理中就没有出现类似的情况,不知道是什么原因呢?2)因为没能成功用包含三条序列的文件跑出来,所以用了包含两条序列的文件,但是我觉得我的碱基种类及数目的统计已经在else:之前的操作中全部完成了,我最后写的else的一串代码并没有发挥作用(最后一行的print没有产生结果),而且我也不是很能理解最后else的作用,望解答,谢谢!这个八十六分钟的视频我居然花了两天才看完的,感觉生信的学习还是任重而道远orz求各路前辈们指教,谢谢!
回复 支持 反对

使用道具 举报

0

主题

3

帖子

39

积分

新手上路

Rank: 1

积分
39
发表于 2017-7-12 01:59:56 | 显示全部楼层
旭日早升 发表于 2017-1-15 17:31
201-python-无名

自己的思考解答:学习了第一题之后,感觉python编程处理生物问题挺有意思的。第二题直接 ...

请问层主是在哪里看的视频,我找不到链接啊.....
回复 支持 反对

使用道具 举报

10

主题

52

帖子

559

积分

版主

Rank: 7Rank: 7Rank: 7

积分
559
QQ
发表于 2017-7-12 15:58:17 | 显示全部楼层
芋头 发表于 2017-7-12 01:59
请问层主是在哪里看的视频,我找不到链接啊.....

视频是要购买的,可以联系版主询问。
回复 支持 反对

使用道具 举报

0

主题

3

帖子

39

积分

新手上路

Rank: 1

积分
39
发表于 2017-7-13 03:41:32 | 显示全部楼层
旭日早升 发表于 2017-7-12 15:58
视频是要购买的,可以联系版主询问。

好的,谢啦
回复 支持 反对

使用道具 举报

0

主题

4

帖子

103

积分

注册会员

Rank: 2

积分
103
发表于 2017-9-1 17:18:41 | 显示全部楼层
#!/usr/bin/perl -w

#hg19每条染色体长度,每条染色体N的含量,GC含量。(fasta文件的探索)
use strict;

open IN,"<$ARGV[0]" or die $!;
print "chr\tA\tT\tGC\tN\tlength\tGC(%)\tN(%)\n";

my $id;
my $seq;
my %seq;
while(<IN>){
        chomp;
        if ($_=~/^>/){
                $id = $_;
        }
        else {
                $seq{$id} .= $_;
        }
}
=cut
        foreach $id (sort keys %seq){
        print "$id\n$seq{$id}\n";
        }
打印哈希
=cut
my $GC;
my $seq_len;
my $GC_rat;
my @GC;
my @A;
my @T;
my @N;
my $A;
my $T;
my $N;

foreach $id (sort keys %seq){
         $seq_len = length $seq{$id};


         @GC = ($seq{$id} =~/(G|C|g|c)/g);
         $GC = $#GC +1;
         @A = ($seq{$id} =~/(A|a)/g);
         $A = $#A +1;
         @T = ($seq{$id} =~/(T|t)/g);
         $T = $#T +1;
         @N = ($seq{$id} =~/(N|n)/g);
         $N = $#N +1;       
     
         my $percent_GC = ($GC)/($A+$T+$GC);
         my $percent_N = ($N)/($A+$T+$GC);
        print "$id\t$A\t$T\t$GC\t$N\t$seq_len\t$percent_GC\t$percent_N\n";
}


close IN;

E:\资料\perl\生信技能树\第二题
回复 支持 反对

使用道具 举报

3

主题

8

帖子

197

积分

注册会员

Rank: 2

积分
197
发表于 2017-12-27 19:30:50 | 显示全部楼层
[mw_shl_code=python,true]        
#! /bin/python
import pysam
import os
os.chdir("/home/liuzh/gpfs/zhangxt/practice")
hg19 = pysam.FastaFile("hg19.fa")
for chr in hg19.references:
        seq = hg19.fetch(chr)
        seqLen = len(seq)
        N = seq.count("N")
        GC = seq.count("G")+seq.count("C")+seq.count("g")+seq.count("c")
        print(chr,"N%","%.2f"%(N/seqLen),"GC%","%.2f"%(GC/(seqLen-N)))
[/mw_shl_code]

回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2020-1-21 06:20 , Processed in 0.024751 second(s), 21 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.