下载数据（人类）

NCBI

由于我需要获取剪切位点两边的序列，那么我需要下载参考基因组数据和注释文件。参考基因组下载常用的有ncbi、ucsc和ensemble。下图是参考基因组版本对应信息。

我是从NCBI下载

链接（https://www.ncbi.nlm.nih.gov/）
点击搜索之后，就可以在页面中找到了。

GenCode

我的注释文件是在GenCode下载，下面为版本信息。跟上面下载的参考基因组版本对应就行。

然后，如下图所示，我现在的该注释文件，GTF文件格式。

下载完之后，我们可以发现注释文件的规律。一个基因号下面是多个转录本，每个转录本下面是外显子，外显子又会分为UTR和CDS。如下图所示：

注释文件中是没有内含子信息，但是通过上图我们可以发现，基因中除了外显子就是内含子。所以我们可以根据注释文件中基因的位置和外显子的位置来确定内含子的位置。

bedtools

获取内含子位置

通过下述代码获取内含子位置，注意修改注释文件名称。

# 安装依赖
$ conda install bedtools
$ conda install bedops


# Gene annotation file gencode.v19.annotation.gtf
GENE_ANNOT_FILE='gencode.v19.annotation.gtf'# Extract genes
awk '$3 == "gene"' $GENE_ANNOT_FILE > genes_temp.gtf# Use BEDops convert2bed
# convert2bed won't work if the "transcript_id" field is not there. This makes an empty "transcript_id" field so convert2bed is happy.
awk '{ if ($0 ~ "transcript_id") print $0; else print $0" transcript_id \"\";"; }' genes_temp.gtf > genes.gtf#Extract exons (exons always have transcript_id field, no need to repeat above)
awk '$3 == "exon"' $GENE_ANNOT_FILE > exons.gtf#Convert both to bed format so we can use bedtools
convert2bed -i gtf < exons.gtf > exons.bed
convert2bed -i gtf < genes.gtf > genes.bed# Remove exons from the genes.bed file, enforce strand specificy (-s)
bedtools subtract -a genes.bed -b exons.bed -s > introns.bed# Convert intron bed file to GTF file
# 最终结果
awk '{print $1"\t"$7"\t""intron""\t"($2+1)"\t"$3"\t"$5"\t"$6"\t"$9"\t"(substr($0, index($0,$10)))}' introns.bed > introns.gtf# combine with original annotations file
cat gencode.v19.annotation.gtf introns.gtf > gencode.v19.introns.gtf# remove generated temporary files
rm genes.gtf exons.gtf exons.bed genes.bed introns.bed

获取内含子序列

利用bedtools 里面的
在获取内含子位置之后，本以为获取内含子序列会非常快速。但是，被窝瞎搞，浪费了一整天时间。一直注视文件和参考基因组匹配不上，报错：

WARNING. chromosome (chromosome Y) was not found in the FASTA file. Skipping.
在报错之后，应该是参考基因组的描述行没匹配上，而我却在参考基因组上修改，由于文件太大，每次修改都很慢，保存也慢。又要上传服务器运行，浪费了很多时间。以后就在一个小文件测试通过再看吧。记住啦~
后面怎么解决这个错误呢？

# 注释文件如下
chr1    12227   12594   ENSG00000223972.4   .   +  HAVANA  intron  .   gene_id "ENSG00000223972.4"; transcript_id "ENSG00000223972.4"; gene_type "pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1"; level 2; havana_gene "OTTHUMG00000000961.2";
chr1    12721   12974   ENSG00000223972.4   .   +  HAVANA  intron  .   gene_id "ENSG00000223972.4"; transcript_id "ENSG00000223972.4"; gene_type "pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1"; level 2; havana_gene "OTTHUMG00000000961.2";
chr1    13052   13220   ENSG00000223972.4   .   +  HAVANA  intron  .   gene_id "ENSG00000223972.4"; transcript_id "ENSG00000223972.4"; gene_type "pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1"; level 2; havana_gene "OTTHUMG00000000961.2";

# 本来参考基因组fna文件的描述行是下面这样
>NC_000001.10 Homo sapiens chromosome 1, GRCh37.p13 Primary Assembly
# 我尝试了  匹配不上第一次
>NC_000001.10 Homo sapiens chromosome 1, GRCh37.p13 Primary Assembly, chr1
# 我尝试了  匹配不上第二次
>chr1， NC_000001.10 Homo sapiens chromosome 1, GRCh37.p13 Primary Assembly
# 我尝试了  匹配不上第三次；这里也把注释文件改为chromosome 1
>chromosome 1
># 最终成功的是
>chr1

也就是说，描述行和bed注释文件完全一样才匹配上了。

bedtools getfasta -fi ncbi_dataset/data/GCF_000001405.25/human_ref_GRCh37.fna -bed introns.bed -fo introns.fasta

获取上下游位置和序列

这里获取序列的位置是以上述获得的

# 修改序列长度
# awk '{print $1"\t"($2-300)"\t"($3+300)}' tmp.bed > introns.bed
# -300是指向上游序列再取300 +300是指向下游序列再取300
# 这里有个问题，当$2-300时候，左边的序列长度居然是301
# 所以这里改为$2-299
awk '{print $1"\t"$7"\t""intron""\t"($2-299)"\t"($3+300)"\t"$5"\t"$6"\t"$9"\t"(substr($0, index($0,$10)))}' introns.bed > example_1.gtf
# 利用convert2bed包将gtf格式转化为bed格式
convert2bed -i gtf < example_1.gtf > example_1.bed
# 获取序列
bedtools getfasta -fi human_ref_GRCh37.fna -bed example_1.bed -fo example_1.fasta

这里虽然获取到了供体和受体对的序列，但是长度可能还会改变。而且还不能作为模型的输入，还有最后一步是提取供体和受体位点侧翼序列拼接在一起（例如文章DeepSplice中长度为120的四个部分序列）

暂时先记录到这里~

参考文章

链接: 下载参考基因组.
链接: 下载注释文件.
链接: UTR和CDS.
链接: bedtools文档.
链接: bedtools使用手册.
链接: bedtools常用命令.
链接: getfasta获取序列.
链接: GTF获取内含子位置.
链接: 获取内含子区域带图.
链接: 获取上下游序列.
链接: bedtools博客教程.
链接: bedtools博客教程2.
链接: DeepSplice.

bedtools从剪切位点两边提取序列相关推荐

基于 PacBio 测序数据的纠错算法评测与剪切位点识别研究
基于 PacBio 测序数据的纠错算法评测与剪切位点识别研究摘要高通量测序技术的产生和发展催生了许多大规模基因测序项目, 如国际千人基因组计划. 英国 UK10K 计划以及中国的百万人群基因组 ...
R语言丨根据VCF文件自动填充对其变异位点并生成序列fa文件
根据VCF文件自动填充对其变异位点并生成序列fa文件首先提出一个问题: 假如有一个重测序结果VCF文件,里面包含了很多个样本在几百个突变位点(snp和iad)的基因型数据,现在想根据这份原始数据,得 ...
根据基因或者蛋白的id提取序列---extract_seq.exe
在面对需要提取某一个基因的序列时,大家会通过什么样的方式提取呢?下载具有序列的全部文件,手动一个一个去搜?刚开始,小编也是这样的,那时候从NCBI上一点一点的搜,一点一点的查,经常一查就是一天,关键是 ...
根据ID从FASTA文件中批量提取序列【Python】
根据ID从FASTA文件中批量提取序列[Python] 生信问题记录我的需求 input: FASTA文件,含六千余个蛋白序列.命名为FA.fasta txt文件,经过interpro注释后,筛选出 ...
提取基因结构信息linux,求助：哪位高手知道如何通过基因编号提取序列
求助:哪位高手知道如何通过基因编号提取序列发布时间:2009-05-24 03:12:53来源:红联作者:huangqp [i=s] 本帖最后由 huangqp 于 2009-5-24 03:18 ...
R语言实现GWAS结果显著SNP位点归类提取与变异类型转化
GWAS结果显著SNP位点归类提取与变异类型转化根据GWAS得到的Rresult文件信息,能够找出每个snp位点对应的显著性情况和基因变异信息,接下来,需要根据表格中的信息进行归纳总结,对不同显著性 ...
linux批量筛选序列变异位点,使用bedtools获取指定坐标上下游的序列
前些天给学徒演示了猪狗的参考基因组构建索引就顺便布置了作业,有意思的是她下载的时候,在两个参考基因组文件里面犹豫不决: : The systematic name of the species. : ...
nodejs剪切视频，提取音频，上传播放
简单说说实现方案首先要有演唱会的链接,使用ibili 这个库下载视频,也可以自己抓取视频链接请求下载,这里有很多方法. 将视频保存在本地后,整理出每一首歌曲对应的时分秒.我找的这个视频在某站评论中已 ...
python文本提取序列信息_从fasta文件中通过头中的ID号提取序列
accessionids.txt是否只包含四位数代码?在如果是,请将accessorID更改为:accessorID = accessorIDWithArrow[1:5] 一些方法可以让这更像Pyt ...

bedtools从剪切位点两边提取序列

目录