搜索
楼主: Jimmy

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

  [复制链接]

0

主题

2

帖子

267

积分

中级会员

Rank: 3Rank: 3

积分
267
发表于 2017-1-15 21:00:40 | 显示全部楼层
本帖最后由 钻石狗 于 2017-1-15 21:02 编辑

089-python——钻石狗
学习了老师教的pysam模块
import pysam
hg19 = pysam.FastaFile("genome.fa")
for seqName in hg19.references:
        seq = hg19.fetch(seqName)
        seqLen = len(seq)
        N = seq.count("N")
        GC = seq.count("G") + seq.count("g") + seq.count("C") + seq.count("c")
        print(seqName,seqLen, (N/seqLen),(GC/seqLen))

因为要放假了事情多,另外一种方法还没有深入学习。以后尽快补上。
回复 支持 反对

使用道具 举报

2

主题

21

帖子

243

积分

中级会员

Rank: 3Rank: 3

积分
243
发表于 2017-1-15 21:37:47 | 显示全部楼层
147-R+perl+bash

本来以为今天能够完成自身工作,但是没有完成,作业明晚好好思考,今天晚上先看看视频,明天晚上自己再好好思考并完成
回复 支持 反对

使用道具 举报

0

主题

3

帖子

24

积分

新手上路

Rank: 1

积分
24
发表于 2017-1-15 21:50:33 | 显示全部楼层
016 Python+R
最近在忙着考试,抽空看了视频,25号之前补上作业。
回复 支持 反对

使用道具 举报

1

主题

6

帖子

434

积分

中级会员

Rank: 3Rank: 3

积分
434
发表于 2017-1-15 21:51:17 | 显示全部楼层
041-python
计算fasta文件中的N和GC含量,看了视频教程思路基本明白了,仿照着代码用测试文件运行了一遍,代码和结果如下:
[Python] 纯文本查看 复制代码
import os
from collections import OrderedDict
os.chdir("C:/Users/de/Desktop/")

chr_dict = OrderedDict()  #定义一个有序字典
temp_dict = OrderedDict()
temp_chr = ""

with open("test.fasta.txt","r") as test: #打开文件
    for line in test:
        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():
    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)))  #定义输出结果保留两位小数

结果如图:
QQ截图20170115214750.jpg

通过本次练习,体会到python的语言魅力,编程基础薄弱,还在学习中,为有趣的事情要坚持下去。
回复 支持 反对

使用道具 举报

0

主题

4

帖子

59

积分

注册会员

Rank: 2

积分
59
发表于 2017-1-16 10:32:15 | 显示全部楼层
本帖最后由 ws34sd 于 2017-1-16 10:55 编辑

090-py+R
[Python] 纯文本查看 复制代码
import re
import sys
import copy
import os

from collections import OrderedDict 

#os.chdir("D:\bio")

input_file_name="sample.txt"

file=open(input_file_name,"r")

temstring=""
chr_dict=OrderedDict()

for line in file:
        line=line.strip()
        if line.startswith(">"):
                temstring=line
                chr_dict[temstring]=""
        else:
                chr_dict[temstring]+=line
#print(chr_dict)

for seqName,seq in chr_dict.items():
        seqlen=len(seq)
        N_ratio=seq.count("N")/seqlen
        GC_number=seq.upper().count("G")+seq.upper().count("C")
        GC_ration=GC_number/(seqlen-seq.count("N"))
        print(seqName,N_ratio,GC_ration)



>chr_1 0.09090909090909091 0.425
>chr_2 0.11764705882352941 0.43333333333333335
>chr_3 0.12121212121212122 0.4482758620689655
>chr_4 0.12 0.4090909090909091



回复 支持 反对

使用道具 举报

3

主题

13

帖子

104

积分

注册会员

Rank: 2

积分
104
发表于 2017-1-16 11:41:36 | 显示全部楼层
本帖最后由 猪兔子2号 于 2017-1-16 15:38 编辑
sxuan 发表于 2017-1-9 21:46
这个upper用后面有啥区别吗

perl :007
[Perl] 纯文本查看 复制代码
[mw_shl_code=perl,true]#!/usr/bin/perl
#
#
while(<>){
if (/^>/){chomp;$chr=$_;$hash{$chr}}else{chomp;
$count{$chr."A"}+=($_=~tr/[Aa]//);
$count{$chr."C"}+=($_=~tr/[Cc]//);
$count{$chr."G"}+=($_=~tr/[Gg]//);
$count{$chr."T"}+=($_=~tr/[Tt]//);
$count{$chr."N"}+=($_=~tr/[Nn]//);
$count{$chr}=$count{$chr."A"}+$count{$chr."T"}+$count{$chr."G"}+$count{$chr."C"};
                                         }
         }
foreach $i(1..22,X,Y,M){
if ($count{">chr".$i}){$sum=$count{">chr".$i};}
if($count{">chr".$i."N"}){$N=$count{">chr".$i."N"}}
if($count{">chr".$i."G"}||$count{">chr".$i."C"}){$GC=$count{">chr".$i."G"}+$count{">chr".$i."C"}};
$ratio=$GC/$sum;
print "chr$i\t$sum\t$N\t$GC\t$ratio\n";

}
[/mw_shl_code]

这里统计了每条染色体的不带N的碱基数,N的碱基数,GC的碱基数以及GC比例。
1 chr1    225280621       23970000        94040974        0.417439252353623
      2 chr2    238204518       4994855 95862507        0.402437820259983
      3 chr3    194797135       3225295 77323307        0.396942732242956
      4 chr4    187661676       3492600 71776628        0.382478881836268
      5 chr5    177695260       3220000 70218569        0.395162870410837
      6 chr6    167395066       3720001 66306710        0.396109106346062
      7 chr7    155353663       3785000 63308649        0.407513075504374
      8 chr8    142888922       3475100 57406604        0.401756855580449
      9 chr9    120143431       21070000        49639471        0.413168415341826
     10 chr10   131314738       4220009 54607071        0.415848760251115
     11 chr11   131129516       3877000 54504836        0.415656502537537
     12 chr12   130481393       3370502 53252045        0.40811983820559
     13 chr13   95589878        19580000        36827474        0.385265414817247
     14 chr14   88289540        19060000        36099079        0.408871526570418
     15 chr15   81694766        20836626        34475969        0.42200952016926
     16 chr16   78884753        11470000        35332028        0.447894259109869
     17 chr17   77795210        3400000 35428296        0.45540459367614
     18 chr18   74657229        3420019 29702356        0.397849697850425
     19 chr19   55808983        3320000 26989400        0.483603150410392
     20 chr20   59505520        3520000 26257240        0.441257214456743
     21 chr21   35106642        13023253        14334933        0.408325381846546
     22 chr22   34894545        16410021        16745219        0.479880709148092
     23 chrX    151100560       4170000 59679184        0.394963354205967
     24 chrY    25653566        33720000        10252459        0.399650442359553
     25 chrM    16571   33720000        7372    0.444873574316577
     26 0.400864728838764


暂时我也沿用了hash的,只是主要格式跟群主不太一样,但是大体的主题原理是一样的,还可以采用split的方法,但是好像没有这个方便。下面统计了gene的GCratio
[Perl] 纯文本查看 复制代码
while(<>){
if (/^>/){chomp;$chr=$_;$hash{$chr}}else{chomp;
$count{$chr."A"}+=($_=~tr/[Aa]//);
$count{$chr."C"}+=($_=~tr/[Cc]//);
$count{$chr."G"}+=($_=~tr/[Gg]//);
$count{$chr."T"}+=($_=~tr/[Tt]//);
$count{$chr."N"}+=($_=~tr/[Nn]//);
$count{$chr}=$count{$chr."A"}+$count{$chr."T"}+$count{$chr."G"}+$count{$chr."C"};
$ratio{$chr}=($count{$chr."C"}+$count{$chr."G"})/$count{$chr};
}
}
$i=0;
foreach $ratio(keys %ratio){

$i++;
$sumratio+=$ratio{$ratio};
print $ratio."\t".$ratio{$ratio}."\n";


}
$averatio=$sumratio/$i;
print $averatio."\n";


所有gene的GCratio是0.494290110527937,上面同得到的genome的GC含量是0.400864728838764,可以看出exon区域的GC含量是高一点。




回复 支持 反对

使用道具 举报

0

主题

22

帖子

183

积分

注册会员

Rank: 2

积分
183
发表于 2017-1-16 18:35:43 | 显示全部楼层
ybb7 发表于 2017-1-14 11:22
004-R-十
心得
1:利用tapply实现类似于excel里的数据透视表的功能tapply(myseq,index,paste,collapse="")

[Tt]表示用正则表达式匹配字母T和t,虽然和你表述的不区分大小写的结果一致,但是意思不同。
[]表示的是对list取其中的元素,因为gregexpr返回的结果是list数据结构,所以用[]表示返回第i个元素的意思。具体可以google关键词“R 字符串 正则”。
参考链接
http://rstudio-pubs-static.s3.am ... b4fbd2ea5ea162.html
http://blog.qiubio.com:8080/archives/2876
http://yphuang.github.io/blog/20 ... gs-processing-in-R/
回复 支持 反对

使用道具 举报

0

主题

22

帖子

183

积分

注册会员

Rank: 2

积分
183
发表于 2017-1-16 18:40:13 | 显示全部楼层
ybb7 发表于 2017-1-14 11:22
004-R-十
心得
1:利用tapply实现类似于excel里的数据透视表的功能tapply(myseq,index,paste,collapse="")

[Tt]的意思是表示用正则表达式匹配字母T和t,虽然结果和你想象的忽略大小写是一样的,但是含义不同。
[]的含义是对list取第i个元素,因为前面gregexpr正则表示返回的结果是list结构。
google搜索关键词“R+字符串+正则”
具体可以参考
http://rstudio-pubs-static.s3.am ... b4fbd2ea5ea162.html
http://blog.qiubio.com:8080/archives/2876
http://yphuang.github.io/blog/20 ... gs-processing-in-R/
回复 支持 反对

使用道具 举报

0

主题

22

帖子

183

积分

注册会员

Rank: 2

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

tr///,在perl里面实际的操作结果应该是字符替换,比如tr/Aa//,就是将Aa字符转换为“”(无字符,结果有点像是删除字符),但是他的返回值是替换字符的数量,这一点要注意。
参考链接
http://odyssey.blog.sohu.com/269968330.html
回复 支持 反对

使用道具 举报

0

主题

22

帖子

183

积分

注册会员

Rank: 2

积分
183
发表于 2017-1-16 19:15:20 | 显示全部楼层
Andyhigh 发表于 2017-1-15 15:01
下载了R和python的视频,想着对R熟悉些,就把R的视频看了下,发现代码等看不清;

于是又来论坛看大家的 ...

R的视频会非常不清晰吗?理应看得清楚,看此贴上其他同学的回复,有贴出同样代码的。同时有位同学自己写的R代码也非常棒,在前面沙发的位置上。
回复 支持 反对

使用道具 举报

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

本版积分规则

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

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

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.