搜索
查看: 2991|回复: 14

生信编程直播第9题:根据指定染色体及坐标得到参考碱基

[复制链接]

633

主题

1172

帖子

3947

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3947
发表于 2017-1-15 21:04:48 | 显示全部楼层 |阅读模式
参考基因组,假设是hg19吧!
指定染色体及坐标,假设是"chr5","8397384"
那么如何写程序得到 这个坐标以及它上下一个碱基呢?
print &get_context("chr5","8397384");
可以看到我写的这个函数,做到了取第五条染色体的8397384位点的上下一个碱基,在UCSC里面也可以验证一下。
当然,要做出批量的!
因为我们是根据vcf文件来做这件事情。
而VCF文件里面记录了所有的变异位点的坐标,我们需要知道上下文来做mutation signature的分析。
我把这个需求拆解开来,希望对你们有帮助!

当然,考虑到很多人的机器hold不住hg19这个大基因组,可以用一个小的fasta文件作为例子:
比如,基因组是:
[AppleScript] 纯文本查看 复制代码
>chr_1
ATCGTCGaaAATGAANccNNttGTA
AGGTCTNAAccAAttGggG
>chr_2
ATCGAATGATCGANNNGccTA
AGGTCTNAAAAGG
>chr_3
ATCGTCGANNNGTAATggGA
AGGTCTNAAAAGG
>chr_4
ATCGTCaaaGANNAATGANGgggTA

指定坐标是 3号染色体的第6个碱基,用程序算,是什么碱基,不允许用肉眼看!

                  http://www.bio-info-trainee.com/?p=1049



上一篇:R plot 之 bubble plot
下一篇:对snp进行注释并格式化输出-高级难度题目
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

10

主题

50

帖子

219

积分

中级会员

Rank: 3Rank: 3

积分
219
QQ
发表于 2017-1-15 22:35:38 | 显示全部楼层
这个真的很实用
苛求远离完美
回复 支持 反对

使用道具 举报

10

主题

50

帖子

498

积分

版主

Rank: 7Rank: 7Rank: 7

积分
498
QQ
发表于 2017-3-20 18:44:04 | 显示全部楼层
本帖最后由 旭日早升 于 2017-3-20 18:45 编辑

201-python-无名
这次作业是提取特定染色体上特定位置的碱基信息,这道题有点类似第二次作业,也是对基因组序列探索,因而我使用了第二题李老师讲解的方法。
解题思路:1、导入pysam包,2、用pysam包读取fasta文件,3、输出对应信息。
我的代码
[Python] 纯文本查看 复制代码
###import module
import os, pysam
###change directory if necessary
os.chdir("./")
###read fasta file use pysam
hg19 = pysam.FastaFile("chrome.fasta") 
dir(hg19) #list methods
###print the 6th on chr_3
print(hg19.fetch("chr_3")[5])

结果展示
C



回复 支持 反对

使用道具 举报

0

主题

3

帖子

138

积分

注册会员

Rank: 2

积分
138
发表于 2017-4-6 21:13:50 | 显示全部楼层
直接samtools faidx hg19.reference chr1:12953463-12953463
hg19.reference 要先建立索引 直接samtools faidx hg19.reference
回复 支持 反对

使用道具 举报

64

主题

138

帖子

681

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
681
发表于 2017-4-7 09:45:37 | 显示全部楼层
R语言实现的,用得是测试数据。
首先读取fasta文件
读取序列
取第3条染色体
然后用字符串截取读取第6个字符。
比较傻
[Python] 纯文本查看 复制代码
library(ShortRead)
library(Biostrings)
setwd('D:\\')#设置文件存放路径
#指定文件名称
file="test.fasta"
#读入fasta文件
seqs<-readFasta(file)#读取fasta文件
x<-sread(seqs)读取序列
x<-list(x)
x[[1]][3]#取第3条染色体
substr(x[[1]][3], 6, 6)#读取第6个字符
回复 支持 反对

使用道具 举报

2

主题

34

帖子

553

积分

高级会员

Rank: 4

积分
553
发表于 2017-4-13 00:38:44 | 显示全部楼层
本帖最后由 x2yline 于 2017-4-13 01:07 编辑

077-x2yline
用python3中buffer的方式读取文件,使用下载并合并好的hg19.fa文件计算
The base of chr3:6 is N
The base of chr5:8397384 is C


代码和第二题的代码类似,
[Python] 纯文本查看 复制代码
file = r'E:\r\biotrainee_demo\class2\hg19.fa'
chromosome = 5#int(input("Please enter the chromosome:\n"))
base_pos = 8397384#int(input("Please enter the base position:\n"))
with open(file, 'r') as f:
    head = ''
    got = 0
    while True:
        base_begin_pos = 0
        buffer = head + f.read(2048*2048)
        if '>' in buffer:
            print(buffer[buffer.find('>')+1:buffer.find('>')+6].strip(), end='\r')
        if buffer.rfind('>') > buffer.rfind('\n'):
            buffer = buffer[:buffer.rfind('>')]
            head = buffer[buffer.rfind('>'):]
        else:
            buffer = buffer
            head = ''
        if ('>chr'+ str(chromosome)+'\n') in buffer:
            got = 1
            chr_begin = buffer[buffer.find('>chr'+ str(chromosome)+'\n')+len('>chr'+ str(chromosome))+1:].replace('\n','').replace(' ','')
            base_begin_pos += len(chr_begin)
            if base_begin_pos >=base_pos:
                end = True
                base_value = chr_begin[base_pos-1]
                break
            else:
                chr_skip1 = f.read(base_pos-base_begin_pos-1).replace('\n','').replace(' ','')
                chr_skip2 = f.read(base_pos-base_begin_pos-1).replace('\n','').replace(' ','')
                base_value = chr_skip2[base_pos-base_begin_pos-len(chr_skip1)-1]
                end = True
                if '>' in chr_skip1 or '>' in chr_skip2[:base_pos-base_begin_pos-len(chr_skip1)]:
                    print("\nThe position do not exist!!!!\n")
                    end = False
                    break
        if got == 1 or len(buffer)==0:
            break
    if end:
        print("\nThe base of chr{}:{} is {}".format(chromosome, base_pos, base_value))
如下:

回复 支持 反对

使用道具 举报

0

主题

5

帖子

99

积分

注册会员

Rank: 2

积分
99
发表于 2017-4-24 20:47:46 | 显示全部楼层
[Python] 纯文本查看 复制代码
import os
from collections import OrderedDict
os.chdir("E:\\work\\biopython\\biotrainee\\lesson9")

gene_seq = OrderedDict()

def get_seq(file,id,pos=0):
    with open(file,'rt') as f:
        for line in f:
            if line.startswith('>'):
                ID = line.strip('\n')[1:]
                gene_seq[ID] = ''
            else:
                gene_seq[ID] += line
        
        if pos:
             seq = gene_seq[id][pos-1]
        else:
             seq = gene_seq[id]
        
        print(file,id,pos,seq)
        
get_seq("test.fa",'Chr1',50)
get_seq("test.fa",'Chr1',52)
回复 支持 反对

使用道具 举报

6

主题

23

帖子

185

积分

注册会员

Rank: 2

积分
185
发表于 2017-5-14 22:39:02 | 显示全部楼层
本帖最后由 qin_qinyang 于 2017-5-14 22:40 编辑

123_python_qin
[AppleScript] 纯文本查看 复制代码
import os
from collections import OrderedDict
os.chdir("/Users/qin/Desktop/biotree/python/lession_9")
sequence = OrderedDict()
with open ("file.test.txt", 'rt') as f:
        for line in f:
                line = line.strip()
                line = line.upper()
                if line.startswith('>'):
                        id = line.split('_')[1]
                        sequence[id]=''
                else:
                        sequence[id]+=line
position = 'chr1:6'
find_chr = position.split(':')[0][3:]
find_pos = position.split(':')[1]
for id, value in sequence.items():
        if find_chr == id:
                print(id,value[int(find_pos)-1])

1 C


回复 支持 反对

使用道具 举报

0

主题

16

帖子

145

积分

注册会员

Rank: 2

积分
145
发表于 2017-5-15 14:32:46 | 显示全部楼层
继续补作业

整个过程不复杂,标准流程读取FASTA文件,做成字典。
通过染色体号码重新组装成key去读取序列,最后取碱基。

[Python] 纯文本查看 复制代码
#!/usr/bin/env python
# encoding=utf-8

def findbase(filename, chrom, pos):
    fi = open(filename, 'r')
    ID = ''
    genedic = {}
    pos = int(pos)
    for line in fi:
        line = line.strip('\n')
        if line.startswith('>'):
            ID = line
            genedic[ID] = ''
        else:
            genedic[ID] += line
    key = '>chr_' + str(chrom)
    base = genedic[key][pos - 1]
    print(base)


findbase('fa1.fa', 3, 6)


输出C
回复 支持 反对

使用道具 举报

0

主题

15

帖子

147

积分

注册会员

Rank: 2

积分
147
发表于 2017-6-1 15:37:21 | 显示全部楼层
毕业的事情,整了这么久,终于过来继续做题了,实在是对不起自己了
代码贴一下,用的比较笨的办法啦

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

use warnings;
use strict;

open A,$ARGV[0] or die $!;
$/ = ">";
<A>;
my %hash;
while(<A>){
	chomp;
	my($info,@seq) = split(/\n/,$_,2);	
	$hash{$info} = join "",@seq;
}
close A;
foreach my $k(sort keys %hash){
	if ($k eq "chr_3"){
		my $seq = $hash{$k};
		$seq =~ s/\n//g;
		print "$seq\n";
		my $goal = substr($seq,5,1);
		print "$goal\n";
	}else{
		next;
	}
}




分别展示了全长和第6个碱基




本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2018-2-23 06:06 , Processed in 0.129360 second(s), 30 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.