搜索
楼主: Jimmy

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

  [复制链接]

0

主题

12

帖子

165

积分

注册会员

Rank: 2

积分
165
发表于 2017-1-16 20:03:24 | 显示全部楼层
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl
#use warnings;
$a=$ARGV[0];
open IN, "$a" or die "can not open";
print "chr\tLength\(bp\)\tA\(bp\)\tT\(bp\)\tC\(bp\)\tG\(bp\)\tN\(bp\)\tN(%)\tGC(%)\n";
while(<IN>){
    chomp;
    if ($_=~/^>(chr\S+)/){
        $chr=$1;
    }
    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 $chr (sort keys %count_A){
    $len=($count_A{$chr} + $count_T{$chr} + $count_C{$chr} + $count_G{$chr} + $count_N{$chr});
    $GC=(($count_C{$chr}+$count_G{$chr})/$len)*100;
    $n=($count_N{$chr}/$len)*100;
    printf "%s\t\t%d\t%d\t%d\t%d\t%d\t%d\t%.2f\t%.2f\n","$chr","$len","$count_A{$chr}","$count_T{$chr}","$count_C{$chr}","$count_G{$chr}","$count_N{$chr}","$n","$GC";
}
close IN

用的测试数据,结果如下:

yzk84965137_1484568126321_12.png
回复 支持 反对

使用道具 举报

0

主题

4

帖子

97

积分

注册会员

Rank: 2

积分
97
发表于 2017-1-16 20:47:11 | 显示全部楼层
学员编号:099-perl+R+shell
最近一直在忙着期末的事情,缺的作业会一点点补上。第二题在之前刚接触perl的时候就被师姐要求写过,整体来说基本没啥大问题,之前也是用的tr进行计数,这次将染色体进行了批量的处理,主要的改变就是将染色体号当做哈希的key来用哈希进行计数,代码如下:
[Perl] 纯文本查看 复制代码
#! usr/bin/perl 
use warnings;use strict;
use Getopt::Long;
my ($input,$output,$id,%A_count,%T_count,%G_count,%C_count,%N_count,$length,$GC);
GetOptions('input=s' => \$input,
			'output=s' => \$output);
open (FILE,"<$input")|| die"$!";
open (OUT,">$output")|| die "$!";
$length=0;
	while (my $line=<FILE>){
		chomp $line;
		if ($line=~/^>chr/){
			$id=$line;
			$id=~s/>chr/chr/;
		}else{
			$A_count{$id}+=($line=~tr/Aa//);
			$T_count{$id}+=($line=~tr/Tt//);
			$G_count{$id}+=($line=~tr/Gg//);
			$C_count{$id}+=($line=~tr/Cc//);
			$N_count{$id}+=($line=~tr/Nn//);
		}
	}
	foreach my $key(sort keys %A_count){
		$length = $A_count{$key}+$T_count{$key}+$G_count{$key}+$C_count{$key}+$N_count{$key};
		$GC = ($G_count{$key}+$C_count{$key})/($A_count{$key}+$T_count{$key}+$G_count{$key}+$C_count{$key})*100;
		print OUT "chrom:$key\nA_count\tT_count\tG_count\tC_count\tN_count\n";
		print OUT "$A_count{$key}\t$T_count{$key}\t$G_count{$key}\t$C_count{$key}\t$N_count{$key}\n";
		print OUT "The $key chrom length(bp):$length\nThe $key chromGC content(%) is $GC\n";
	}
close (FILE);
close (OUT);


结果图

结果图
回复 支持 反对

使用道具 举报

2

主题

21

帖子

243

积分

中级会员

Rank: 3Rank: 3

积分
243
发表于 2017-1-16 21:25:45 | 显示全部楼层
本帖最后由 AdaWon 于 2017-1-16 21:27 编辑

147-R+perl+bash


###待解决番茄的每条染色体长度,每条染色体A/T/G/C/N的含量,GC含量。
###GC含量是在DNA4种碱基中,鸟嘌呤(G)和胞嘧啶(C)所占的比率称为GC含量。
###以下是部分fasta序列
>SL2.50ch00
AATAATAATAATAATAATAATAATAAATAAATAAATAAATAATAATAATAATAATAATAATAAATAAATA
AATAAATAAATAAATAAATAAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAA
TAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAAA
AATAATAATAATAATAATAATAATAATAATAATCATCATCATCATCATCATCATCATCATCATCATCATA
ATAATAATAATCATAATAAGGATAGATGATCTTTTCAACTTGCTTCATGGTATGTCCTTTGGTTTGAAAA
ATTCACAAGCAGCATTCATGGATCTTATGAATAGAGTCTTCAAACCTTATTTAGATCTATTTCTCTTCGT
CTTCATTGATGACATACTGGTTTATTCAAGAACTTAAGAAGATAATGCTCGTCACTTCAGGATTGTTCTG
CAGACTTCAAATGATAAAGAGTTATATTCCAAGTTCTCCAAGTTGGAGTTCGGCTTAAGTTCGTGGCATT
CTCGGGCCACATTGTTTCTGGTGATGGGATTACGTTGATACCCAAAAGATAGAAGCATTTCAGAGTTGGC
CTAGACCCACATCTCCAATAGATATAAAGAGTTTCTTGGGTTAAGTTGGCTATTATGGGAGGTTTGTTGA

###单行命令的收集
[Perl] 纯文本查看 复制代码
perl -alne '{if(/^>/){$chr=$_}else{ $A_count{$chr}+=($_=~tr/Aa//); $T_count{$chr}+=($_=~tr/Tt//);$C_count{$chr}+=($_=~tr/Cc//); $G_count{$chr}+=($_=~tr/Gg//); $N_count{$chr}+=($_=~tr/Nn//); }}END{print "$_\t$A_count{$_}\t$T_count{$_}\t$C_count{$_}\t$G_count{$_}\t$N_count{$_}" foreach sort keys  %N_count}' test.fa


[Perl] 纯文本查看 复制代码
#!/usr/bin/perl

open IN,"<$ARGV[0]";
$/ = ">";<IN>;  #这一行表示什么含义?
print "chr\tN(%)\tGC(%)\n";
while(<IN>){
        s/>//;
        my($chr,$seq) = split; #split函数把字符串进行分割并把分割后的结果放入数组中,这里split后面未加任何分隔符?不理解
        $len = length $seq;
        $n = ($seq =~ s/N/N/ig); #全局替换,返回的是次数!考虑大小写
        $g = ($seq =~ s/G/G/ig);
        $c = ($seq =~ s/C/C/ig);
        $N = ($n/$len) * 100;
        $GC = (($g+$c)/$len) * 100;
        printf "%s\t%.2f\t%.2f\n","$chr","$N","$GC";
}

结果如下:
chr        N(%)        GC(%)
SL2.50ch00        0.00        0.00
SL2.50ch01        0.00        50.00
SL2.50ch02        0.00        37.14
SL2.50ch03        0.00        35.71
SL2.50ch04        0.00        52.86
SL2.50ch05        0.00        48.57
SL2.50ch06        0.00        32.86
SL2.50ch07        0.00        42.86
SL2.50ch08        0.00        28.57
SL2.50ch09        0.00        60.00
SL2.50ch10        0.00        38.57
SL2.50ch11        0.00        48.57
SL2.50ch12        0.00        40.00



为什么结果这么奇怪?而且计算速度极快,但是不太明白这段代码的部分地方


[Perl] 纯文本查看 复制代码
#!/usr/bin/perl
print "chr\tA\tT\tC\tG\tN\tLEN\tGC(%)\n";
while (<>){
chomp;
 
 if(/^>/){
 $chr=$_
 }
 else{ 
 $A_count{$chr}+=($_=~tr/Aa//);
 $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{$_}) * 100;
print "$_\t$A_count{$_}\t$T_count{$_}\t$C_count{$_}\t$G_count{$_}\t$N_count{$_}\t$length\t$gc\n" 
 
}


结果如下:
chr        A        T        C        G        N        LEN        GC(%)
>SL2.50ch00        5828429        5781778        3518327        3538187        3139100        21805821        37.802643538734
>SL2.50ch01         28543252        28529725        14550787        14496299        12423381        98543444        33.7285935334255
>SL2.50ch02         15671870        15697606        7928334        7959202        8083432        55340444        33.6194256209005
>SL2.50ch03         20074769        20121498        10320842        10344705        9925850        70787664        33.954865361062
>SL2.50ch04         20046332        20037897        10173967        10190827        6021919        66470942        33.689202884222
>SL2.50ch05         20186419        20239450        10402824        10411937        4634458        65875088        33.9884828095335
>SL2.50ch06         14372797        14408302        7372715        7419137        6178685        49751636        33.9473266339018
>SL2.50ch07         20439152        20390614        10584159        10546887        6084209        68045021        34.1038881156044
>SL2.50ch08         19732817        19677673        10176256        10197942        6081969        65866657        34.0792913396153
>SL2.50ch09         21382831        21364806        11048464        11071682        7614308        72482091        34.1003576459519
>SL2.50ch10         20073867        20078738        10311029        10327550        4736321        65527505        33.949953993329
>SL2.50ch11         16552154        16578958        8592589        8533584        6045240        56302525        34.0769960016742
>SL2.50ch12         20346043        20296909        10585944        10577486        5338821        67145203        34.2414962907876


###上周看完两个半小时的PERL,也在看小骆驼
###R+shell 明晚再战









回复 支持 反对

使用道具 举报

2

主题

21

帖子

243

积分

中级会员

Rank: 3Rank: 3

积分
243
发表于 2017-1-16 21:32:23 | 显示全部楼层
disheng 发表于 2017-1-5 21:25
perl:043
[mw_shl_code=perl,true]#!/usr/bin/perl

请教一下,这两行代码什么意思?
$/ = ">";<IN>;
my($chr,$seq) = split;
回复 支持 反对

使用道具 举报

0

主题

7

帖子

136

积分

注册会员

Rank: 2

积分
136
发表于 2017-1-16 22:19:41 | 显示全部楼层
056-shell-yy

思路:先将fasta文件分割成基本的矩阵形式,行(record)和列(field),然后再逐一对每行的sequence进行总量或者GC含量的统计。
小记:会简单的awk等命令,但基本上没有用过各种变量,正在看shell脚本方面的书,发现还没有看到这几个变量名的解释。所以跟着老师认真学习了这个视频,开始慢慢体会。

数据:用的仍然是jimmy老师给的test.fasta,命名为hg19.fasta
[Shell] 纯文本查看 复制代码
awk 'BEGIN{RS=">";FS="\n"} NR>1{seq="";for (i=2;i<=NF;i++) seq=seq""$i;print $1":"length(seq)}' hg19.fasta 

chr_1:44
chr_2:34
chr_3:33
chr_4:25
刚开始跟着视频中老师的代码来计算自己的数据,结果发现总是提示语法错误,本以为所有都跟老师的代码一样,按理说不应该出现错误的,一个一个代码的检查,结果还是出现了三次报错以后,分别纠正了三次,代码才正确。

总结:每次写代码时应该认真检查语法是否正确,第二还是要多练习代码,很多时候看懂了也以为自己会了,真正到写代码跑程序的时候才发现一堆问题。

回复 支持 反对

使用道具 举报

0

主题

6

帖子

89

积分

注册会员

Rank: 2

积分
89
发表于 2017-1-17 03:07:06 | 显示全部楼层
本帖最后由 banapi 于 2017-2-15 03:32 编辑

158-shell-python-perl-瓜皮
shell
awk命令:经常用于处理大表格,其他方面处理关联数据也经常使用。awk是一种行处理程序,执行awk时,它依次对输入文件中的每一行执行花括号中的代码。
awk简单的使用介绍http://blog.sina.com.cn/s/blog_7b9ace5301014q8o.html

substr函数http://www.cnblogs.com/sunada2005/p/3493941.html

不太会shell,除了视频中老师讲解的代码还模仿了帖子里anlan的代码,使用测试数据进行长度统计和N、CG含量的计算:
[Shell] 纯文本查看 复制代码
awk 'BEGIN{RS=">";FS="\n"}NR>1{seq="";for(i=2;i<=NF;i++)seq=seq""$i;len=length(seq);count_N=0;count_CG=0;for(i=0;i<len;i++){x=substr(seq,i+1,1);if(x~/[nN]/)count_N++;else if(x~/[CcGg]/)count_CG++}print $1":"len" "count_N/len" "count_CG/len}' test

结果
chr_1:44 0.0909091 0.386364
chr_2:34 0.117647 0.382353
chr_3:33 0.121212 0.393939
chr_4:25 0.12 0.36
perl
while(<>) 逐行读取文本
chomp;              #去掉结尾的换行符

$_是perl中的默认变量,读取内容如果不指定的话会放入默认变量中。
tr/// 是对匹配字符进行替换或去冗余等操作,根据不同的参数有不同的效果,当不加参数时是进行匹配字符的计数;
$count=$temp=~tr/A//; #表示计算$temp中出现A的次数,$temp并不改变值

创建了以染色体编号为键,以匹配个数为值的哈希表,所以共有A,T,C,G,N这5个哈希表,之后遍历的时候print的时候也很方便,用foreach读取一个哈希表中所有的键即可。
若是以染色体编号为哈希表,ATCGN为键,会构建多个哈希表,虽然构建的时候代码也是5行,这里是差不多的,但是之后打印不方便。


while(<>){
chomp;
if(/^>/){
$chr = $_;
}
else{
$A_count{$chr} += ($_ =~ tr/Aa//);
$T_count{$chr} += ($_ =~ tr/Tt//);
$C_count{$chr} += ($_ =~ tr/Cc//);
$G_count{$chr} += ($_ =~ tr/Gg//);
$N_count{$chr} += ($_ =~ tr/Nn//);
}
}
foreach(sort keys %A_count){
$length=$A_count{$_}+$T_count{$_}+$C_count{$_}+$G_count{$_}+$N_count{$_};
$CG=($C_count{$_}+$G_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$CG\n";
}

结果:
>chr_1        13        10        7        10        4        0.425
>chr_2        11        6        5        8        4        0.433333333333333
>chr_3        10        6        3        10        4        0.448275862068966
>chr_4        9        4        2        7        3        0.409090909090909


cg含量和上面不容是有没有算N的区别
python
最后进行了浮点数转换,因为是2.7版本,整取相除取整,所以必须先转换浮点数才能相除,否则算出来是0
[Python] 纯文本查看 复制代码
chr_dict = {}
temp = ""

with open("test") as hg19:
        for line in hg19:
                line = line.strip()
                if line.startswith(">"):
                        temp = line
                        chr_dict[temp] = ""
                else:
                        chr_dict[temp] += line
for seqname, seq in chr_dict.items():
        length = len(seq)
        N = seq.count("N")
        CG =seq.count("G") + seq.count("C")+seq.count("c") + seq.count("g")
        print seqname, seq, length, float(CG)/(length-N)

结果:
>chr_4 ATCGTCaaaGANNAATGANGgggTA 25 0.409090909091
>chr_1 ATCGTCGaaAATGAANccNNttGTAAGGTCTNAAccAAttGggG 44 0.425
>chr_2 ATCGAATGATCGANNNGccTAAGGTCTNAAAAGG 34 0.433333333333
>chr_3 ATCGTCGANNNGTAATggGAAGGTCTNAAAAGG 33 0.448275862069









回复 支持 反对

使用道具 举报

633

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
 楼主| 发表于 2017-1-17 09:04:47 | 显示全部楼层
scyxr 发表于 2017-1-16 22:19
056-shell-yy

思路:先将fasta文件分割成基本的矩阵形式,行(record)和列(field),然后再逐一对每行的 ...

刚开始跟着视频中老师的代码来计算自己的数据,结果发现总是提示语法错误,本以为所有都跟老师的代码一样,按理说不应该出现错误的,一个一个代码的检查,结果还是出现了三次报错以后,分别纠正了三次,代码才正确。

总结:每次写代码时应该认真检查语法是否正确,第二还是要多练习代码,很多时候看懂了也以为自己会了,真正到写代码跑程序的时候才发现一堆问题。

这个总结很赞
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

0

主题

12

帖子

165

积分

注册会员

Rank: 2

积分
165
发表于 2017-1-17 09:53:45 | 显示全部楼层
Dorom 发表于 2017-1-16 20:03
[mw_shl_code=perl,true]#!/usr/bin/perl
#use warnings;
$a=$ARGV[0];

学员107-perl+R-Dorom,忘记加了
回复 支持 反对

使用道具 举报

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

[Python] 纯文本查看 复制代码
#第一种方法
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))
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-12-8 14:42 , Processed in 0.036106 second(s), 24 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.