生信编程​9.bedtools从bam文件中提取目的序列

有一些五六年前的学生们都成长为了各个生物信息学相关公司的小领导,而且他们都有了自己的公众号,知乎号,也算是一番人物。最近他们跟我反馈面试找不到或者说很难直接考核筛选到认真干活的生信工程师,挺有意思的。让我想起来了早在生信技能树论坛创立之初我为了引流,而规划的200个生信工程师面试题。值得继续分享:

200个生信工程师面试考题

使用bedtools进行坐标映射:

  1. 将bam文件转换为bed文件

bedtools bamtobed -i ../CA_N_805_RNA_recal.bam > CA_N_805_RNA_recal.bed

结果文件如下

chr1    13146   13296   ST-E00192:663:HL5C3CCXY:4:2221:19055:55561/1    60      +
chr1    13244   13394   ST-E00192:663:HL5C3CCXY:4:2221:19055:55561/2    60      -
chr1    14184   14334   ST-E00192:663:HL5C3CCXY:4:1207:15168:3436/1     60      +

  1. 取出bed文件中与编辑位点有交集的记录,并同时输出两个文件的原始记录

bedtools intersect -a CA_N_805_RNA_recal.bed -b ../CA_N_805_RNA.AtoI_nonAlu_rmHomopolymer.avinput -wo > CA_N_805_RNA.AtoI_nonAlu_rmHomopolymer.bed

结果文件如下

chr1    14462   14612   ST-E00192:663:HL5C3CCXY:4:2205:20811:67709/1    60      +       chr1    14610 14610    T       C       het     495.64  51      0
chr1    14462   14612   ST-E00192:663:HL5C3CCXY:4:2205:20882:67867/1    60      +       chr1    14610 14610    T       C       het     495.64  51      0
chr1    14466   14616   ST-E00192:663:HL5C3CCXY:4:2101:13514:7778/1     60      +       chr1    14610 14610    T       C       het     495.64  51      0

  1. 将bam文件转换为fastq文件

bedtools bamtofastq -i ../CA_N_805_RNA_recal.bam -fq CA_N_805_RNA_recal_1.fq -fq2 CA_N_805_RNA_recal_2.fq

结果文件如下

@ST-E00192:663:HL5C3CCXY:4:2221:19055:55561/1
CTGAGGAAGGAGAAGGGGATGCACTGTTGGGGAGGCAGCTGTAACTCAAAGCCTTAGCCTCTGTTCCCACGAAGGCAGGGCCATCAGGCACCAAAGGGATTCTGCCAGCATAGTGCTCCTGGACCAGTGATACACCCGGCACCCTGTCCT
+
RR7;8;EEBDFCFFCDDDFFBCFCFBCGBDDDFCDCFCCFBCDFCFDEFFCC@FGDCC@FDFABFD@@FC;FFCCCFCDDC?EECFCCBEB?EEEBCCEEFCEAB?EBBEECBBABEC?EADEB?EBBAEECBEB??;CBEB??EABDNT

  1. 将fastq文件转换为fasta文件

perl -alne '{print $_ if $.%4==1|$.%4==2}' CA_N_805_RNA_recal_1.fq | perl -alne '{s/@/>/g;print}' > CA_N_805_RNA_recal_1.fa

perl -alne '{print $_ if $.%4==1|$.%4==2}' CA_N_805_RNA_recal_2.fq | perl -alne '{s/@/>/g;print}' > CA_N_805_RNA_recal_2.fa

结果文件如下

>ST-E00192:663:HL5C3CCXY:4:2221:19055:55561/1
CTGAGGAAGGAGAAGGGGATGCACTGTTGGGGAGGCAGCTGTAACTCAAAGCCTTAGCCTCTGTTCCCACGAAGGCAGGGCCATCAGGCACCAAAGGGATTCTGCCAGCATAGTGCTCCTGGACCAGTGATACACCCGGCACCCTGTCCT
>ST-E00192:663:HL5C3CCXY:4:1207:15168:3436/1
TGTCACTGAAACCTTTTTTGTGGGAGACTATTCCTCCCATCTGCAACAGCTGCCCCTGCTGACTGCCCTTCTCTCCTCCCTCTCATCCCAGAGAAACAGGTCAGCTGGGAGCTTCTGCCCCCACTGCCTAGGGACCAACAGGGGCAGGAG

  1. 合并两个fasta文件

cat CA_N_805_RNA_recal_1.fa CA_N_805_RNA_recal_2.fa > CA_N_805_RNA_recal.fa

  1. 使用python脚本根据bed文件中与编辑位点有交集的记录从fasta文件中提取出序列

import argparse
import pysam

parser = argparse.ArgumentParser()
parser.add_argument('bed',help='input bed file',type=str)
parser.add_argument('fasta',help='input fasta',type=str)

args = parser.parse_args()
bed = args.bed
fasta = args.fasta

fa = pysam.FastaFile(fasta)

with open(bed,'r') as f:
     for i in f.readlines():
         if i.strip() != '':
            arr = i.strip().split('\t')
            if arr[3] in fa.references:
               seq = fa.fetch(arr[3])
               print(">%s \n %s" % (arr[3],seq))

python extract_seq.py CA_N_805_RNA.AtoI_nonAlu_rmHomopolymer.bed CA_N_805_RNA_recal.fa > CA_N_805_RNA_recal_to_blat.fa

得到的目的序列就是之前想通过python脚本输出的序列,用来进行blat比对

文末友情推荐

(0)

相关推荐

  • 6 GATK4完整流程

    0定义变量 source activate wes #GATK=~/biosoft/gatk/gatk-4.1.2.0/gatk ref=/mnt/f/kelly/bioTree/server/wes ...

  • 转录组学习五(reads比对)

    转录组学习一(软件安装)转录组学习二(数据下载)转录组学习三(数据质控)转录组学习四(参考基因组及gtf注释探究)转录组学习五(reads的比对与samtools排序)转录组学习六(reads计数与标 ...

  • 8 比对及找变异步骤的质控

    使用qualimap对wes的比对bam人家总结测序深度和覆盖度ls -lh *raw.vcf-rwxrwxrwx 1 root root 184M Jun 7 10:58 SRR7696207_ra ...

  • 生信编程直播第11题:把文件内容按照染色体分开写出

    问题描述这个需求很常见,因为一般生物信息学数据比较大,比如sam,vcf,或者gtf,bed都是把所有染色体综合在一起的文件.如果想根据染色体把大文件拆分成小的文件呢?比如:ftp://ftp.ncb ...

  • 生信编程5. 根据GTF文件绘制多个转录本结构

    有一些五六年前的学生们都成长为了各个生物信息学相关公司的小领导,而且他们都有了自己的公众号,知乎号,也算是一番人物.最近他们跟我反馈面试找不到或者说很难直接考核筛选到认真干活的生信工程师,挺有意思的. ...

  • 生信编程11.根据gtf格式的基因注释文件提取人所有基因的染色体坐标

    有一些五六年前的学生们都成长为了各个生物信息学相关公司的小领导,而且他们都有了自己的公众号,知乎号,也算是一番人物.最近他们跟我反馈面试找不到或者说很难直接考核筛选到认真干活的生信工程师,挺有意思的. ...

  • 生信编程18.把文件内容按照染色体分开写出

    有一些五六年前的学生们都成长为了各个生物信息学相关公司的小领导,而且他们都有了自己的公众号,知乎号,也算是一番人物.最近他们跟我反馈面试找不到或者说很难直接考核筛选到认真干活的生信工程师,挺有意思的. ...

  • 生信编程直播第七题:写超几何分布检验!

    下载数据 切换到工作目录:cd d/生信技能树-视频直播/第七讲 kegg2gene(第六讲kegg数据解析结果) 暂时不用新的kegg注释数据为了能够统一答案 差异基因list和背景基因list 关 ...

  • 生信编程直播课程优秀学员作业展示1

    题目 人类基因组外显子区域长度 学员:x2yline 具体题目详情请参考生信技能树论坛 题目数据来源为:ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_huma ...

  • 生信编程直播课程优秀学员学习心得及作业展示3

    学习感悟 首先说明一下,我不算是完全从0开始学习,因为生物的知识和python的语言之前都知道一点,但说实话,我的python距离真正的实践还差的很远,也没有常用所以基本忘完. 真的很感谢群主和老师们 ...

  • 生信编程直播课程优秀学员作业展示2

    题目:hg19基因组序列的一些探究 学员:x2yline 具体题目详情请参考生信技能树论坛 数据来源:http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bi ...

  • 生信编程直播第0题-生信编程很简单!

    貌似顺序有点怪异,我们都把1~8题给讲解完了,现在才开始第0题,很明显,这个是临时加入的,因为有很多没有编程基础的朋友也开始跟着我们学习了. 还不知道怎么回事的先查看历史题目: 生物信息学技能面试题( ...