bed文件坐标的调整（awk的用法）

最新推荐文章于 2024-05-31 15:07:04 发布

qq_39306047

最新推荐文章于 2024-05-31 15:07:04 发布

阅读量2k

点赞数 1

分类专栏： Linux

本文链接：https://blog.csdn.net/qq_39306047/article/details/107620022

版权

Linux 专栏收录该内容

78 篇文章 3 订阅

订阅专栏

bed文件坐标的调整（awk的用法）

。annotatePeak的输入也可以是GRanges对象，你如果用R做peak calling的话,直接就可以衔接上ChIPseeker了。

require(ChIPseeker)
f = getSampleFiles()[[4]]
这里我们需要的是一个TxDb对象，这个TxDb就包含了我们需要的各种信息，ChIPseeker会把信息抽取出来，用于注释时使用。

require(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb = TxDb.Hsapiens.UCSC.hg19.knownGene
x = annotatePeak(f, tssRegion=c(-1000, 1000), TxDb=txdb)

Bioconductor里有很多包是针对GRanges对象的，这样方便你在R里做后续的处理，如果你说你不懂这些，只想输出个Excel表格。那么也很容易，用as.data.frame就可以转成data.frame，然后你就可以用write.table输出表格了。

第三种注释
两种注释有时候还不够，我想看peak上下游某个范围内（比如说-5k到5k的距离）都有什么基因，annotatePeak也可以做到。
你只要传个参数说你要这个信息，还有什么区间内，就可以了。

x = annotatePeak(f[[4]], tssRegion=c(-1000, 1000), TxDb=txdb, addFlankGeneInfo=TRUE, flankDistance=5000)输出中多三列: flank_txIds, flank_geneIds和flank_gene_distances，在指定范围内所有的基因都被列出。

awk ‘{print $1+1}’ su [【注意进行加减运算要看，有时候最后有空行也被运算！！】

删除空行
下列 5 种方法可以做到。

◈ sed：过滤和替换文本的流编辑器。
◈ grep：输出匹配到的行。
◈ cat：合并文件并打印内容到标准输出。
◈ tr：替换或删除字符。
◈ awk：awk 工具用于执行 awk 语言编写的程序，专门用于文本处理。
◈ perl：Perl 是一种用于处理文本的编程语言。

$ cat 2daygeek.txt
2daygeek.com is a best Linux blog to learn Linux.

It’s FIVE years old blog.

This website is maintained by Magesh M, it’s licensed under CC BY-NC 4.0.

He got two GIRL babys.

Her names are Tanisha & Renusha.

$ sed ‘/^$/d’ 2daygeek.txt
2daygeek.com is a best Linux blog to learn Linux.
It’s FIVE years old blog.
This website is maintained by Magesh M, it’s licensed under CC BY-NC 4.0.
He got two GIRL babes.
Her names are Tanisha & Renusha.

$ grep . 2daygeek.txt
or
$ grep -Ev “^$” 2daygeek.txt
or
$ grep -v -e ‘^$’ 2daygeek.txt
2daygeek.com is a best Linux blog to learn Linux.
It’s FIVE years old blog.
This website is maintained by Magesh M, it’s licensed under CC BY-NC 4.0.
He got two GIRL babes.
Her names are Tanisha & Renusha.
以下是命令展开的细节：

◈ grep：该命令本身。
◈ .：替换任意字符。
◈ ^：匹配字符串开头。
◈ $：匹配字符串结尾。
◈ E：使用扩展正则匹配模式。
◈ e：使用常规正则匹配模式。
◈ v：反向匹配。
◈ 2daygeek.txt：源文件名。

◈ awk：该命令本身。
◈ //：标记匹配范围。
◈ ^：匹配字符串开头。
◈ $：匹配字符串结尾。
◈ .：匹配任意字符。
◈ !：删除匹配的字符串。
◈ 2daygeek.txt：源文件名。

4）使用awk命令进行简单的数学运算

[linuxmi@linux:~/linuxmi迷]$ awk ‘BEGIN{ a=9;b=8;print"(a + b)=",(a + b)}’
(a + b)= 17
[linuxmi@linux:~/linuxmi迷]$ awk ‘BEGIN{ a=100;b=78;print"(a - b)=",(a - b)}’
(a - b)= 22
[linuxmi@linux:~/linuxmi迷]$ awk ‘BEGIN{ a=80;b=5;print"(a / b)=",(a / b)}’
(a / b)= 16
[linuxmi@linux:~/linuxmi迷]$ awk ‘BEGIN{ a=80;b=5;print"(a ÷ b)=",(a / b)}’
(a ÷ b)= 16
[linuxmi@linux:~/linuxmi迷]$ awk ‘BEGIN{ a=10;b=9;print"(a × b)=",(a * b)}’
(a × b)= 90

https://mp.weixin.qq.com/s?src=11&timestamp=1595850298&ver=2486&signature=pS42jRFE8HChSmC4UkdIHZxr1edkIFN5CRcTFpwv5kNKFucOwlN28v10eKogWe45drNOvMcD0swXPTQJGKHrr2cN3kKl-fvWTyb9EiNmBysA6KU8YdZ4fiY7zk4fSOEQ&new=1数据详见：生信练习题：调整GFF文件中的坐标位置——by徐洲更
在这里插入图片描述

因此：

如果是#开头的注释部分，不处理
如果非#开头的部分，则需要将chr8后边的起始位置数值提取出来，实现将该值与后边两个位置列进行相加-1。具体对于第一行而言，就是使得第4列的30280+25234310-1，第5列的30951+25234310-1。并且在数据处理之后，去掉chr8后边这一串:25234310-25266151.
对于需要处理的数据，思路也比较简单，先按照冒号(😃、短横杠(-)和tab分隔符(\t)将数据分列。分列之后，我们需要加上的数值在第2列，需要被加的值则在第6列和第7列。此外，需要排除第2列和第3列进行输出。

我们一步步来看。

首先，我们先将注释文件挑出来不处理，直接打印：
awk ‘{if(/^#/){print $0}}’ target.gff|head

$ awk ‘{if(/^#/){print $0}}’ target.gff|more

This output was generated with AUGUSTUS (version 3.2.3).

简单解释一下：

//：表示匹配
^#：表示开头为#的
因此/^#/表示匹配开头为#号
$0表示整行
if(/^#/){print $0}串起来就是: 如果是以#开头，那么打印整行。

好，我们继续。

对于不是#开头的行，我们需要对其进行分列，并对6、7列与第2列相加后减1。

awk ‘BEGIN {FS="[-:\t]";OFS="\t"}{if(/^#/){print $0}else{$6=$6+$2-1;$7=$7+$2-1;print $0}}’ target.gff

可以看到，分列成功并进行了相关的数值计算。

This output was generated with AUGUSTUS (version 3.2.3).

AUGUSTUS is a gene prediction tool written by M. Stanke (mario.stanke@uni-greifswald.de),

O. Keller, S. König, L. Gerischer and L. Romoth.

Please cite: Mario Stanke, Mark Diekhans, Robert Baertsch, David Haussler (2008),

Using native and syntenically mapped cDNA alignments to improve de novo gene finding

Bioinformatics 24: 637-644, doi 10.1093/bioinformatics/btn013

No extrinsic information on sequences given.

arabidopsis version. Using default transition matrix.

We have hints for 0 sequences and for 0 of the sequences in the input set.

----- prediction on sequence number 1 (length = 31842, name = chr8:25234310-25266151) -----

Constraints/Hints:

(none)

Predicted genes for sequence number 1 on both strands

start gene g1

chr8 25234310 25266151 AUGUSTUS gene 25237205 25239574 0.03 + . g1
chr8 25234310 25266151 AUGUSTUS transcript 25237205 25239574 0.03 + . g1.t1
chr8 25234310 25266151 AUGUSTUS tss 25237205 25237205 . + . transcript_id “g1.t1”; gene_id “g1”;
chr8 25234310 25266151 AUGUSTUS exon 25237205 25237684 . + . transcript_id “g1.t1”; gene_id “g1”;
chr8 25234310 25266151 AUGUSTUS start_codon 25237500 25237502 . + 0 transcript_id “g1.t1”; gene_id “g1”;
chr8 25234310 25266151 AUGUSTUS intron 25237685 25237767 0.62 + . transcript_id “g1.t1”; gene_id “g1”;

FS是输入字段分隔符，中括号里边的分隔符的关系是“或”。即遇到-或者:或者\t，都作为分隔符分开。
OFS是输出字段分隔符。这里指定的是tab分隔。
BEGIN是为awk提供启动的操作。具体在本案例而言，就是告诉awk，我后边的FS是-:\t，OFS是"\t"。
else后边一串就是将第6列重新赋值为 $2 - 1 ，第 7 列重新赋值为$ 2-1。然后输出整行。
下一个任务就是将第2、3列去掉。

awk ‘BEGIN {FS="[-:\t]";OFS="\t"}{if(/^#/){print $0}else{$6=$6+$2-1;$7=$7+$2-1;$2=$3="";print $0}}’ target.gff|more

This output was generated with AUGUSTUS (version 3.2.3).

AUGUSTUS is a gene prediction tool written by M. Stanke (mario.stanke@uni-greifswald.de),

O. Keller, S. König, L. Gerischer and L. Romoth.

Please cite: Mario Stanke, Mark Diekhans, Robert Baertsch, David Haussler (2008),

Using native and syntenically mapped cDNA alignments to improve de novo gene finding

Bioinformatics 24: 637-644, doi 10.1093/bioinformatics/btn013

No extrinsic information on sequences given.

arabidopsis version. Using default transition matrix.

We have hints for 0 sequences and for 0 of the sequences in the input set.

----- prediction on sequence number 1 (length = 31842, name = chr8:25234310-25266151) -----

Constraints/Hints:

(none)

Predicted genes for sequence number 1 on both strands

start gene g1

chr8 AUGUSTUS gene 25237205 25239574 0.03 + . g1
chr8 AUGUSTUS transcript 25237205 25239574 0.03 + . g1.t1
chr8 AUGUSTUS tss 25237205 25237205 . + . transcript_id “g1.t1”; gene_id “g1”;
chr8 AUGUSTUS exon 25237205 25237684 . + . transcript_id “g1.t1”; gene_id “g1”;

看起来似乎第2、3列似乎是没了。但是其实它们只是为空，我们看不到而已。

我们来仔细看看：

$ awk ‘BEGIN {FS="[-:\t]";OFS="\t"}{if(/^#/){print $0}else{$6=$6+$2-1;$7=$7+$2-1;$2=$3="";print $0}}’ target.gff|awk '!/^#/{print $KaTeX parse error: Expected 'EOF', got '}' at position 2: 0}̲'|head|cat -A c\dots$
chr8^II^IAUGUSTUSItranscript^I25237205I25239574^I0.03I+^I.Ig1.t1$
chr8^II^IAUGUSTUSItss^I25237205I25237205^I.I+^I.Itranscript_id “g1.t1”; gene_id “g1”;$
chr8^II^IAUGUSTUSIexon^I25237205I25237684^I.I+^I.Itranscript_id “g1.t1”; gene_id “g1”;$

为了避免注释行对我们的干扰。我这里输出的时候用awk '!/^#/{print $0}'只输出非注释行。

可以看到，chr8后边有3个^I，这个其实就是tab分隔符的真身。

为了结果更为令人愉悦，我还是决定把这个小刺拔掉。

$ awk ‘BEGIN {FS="[-:\t]";OFS="\t"}{if(/^#/){print $0}else{$6=$6+$2-1;$7=$7+$2-1;$2=$3="";sub(/\t+/,"\t");print}}’ target.gff|awk '!/^#/{print $KaTeX parse error: Expected 'EOF', got '}' at position 2: 0}̲'|head|cat -A c\dots$
chr8^IAUGUSTUSItranscript^I25237205I25239574^I0.03I+^I.Ig1.t1$
chr8^IAUGUSTUSItss^I25237205I25237205^I.I+^I.Itranscript_id “g1.t1”; gene_id “g1”;$
chr8^IAUGUSTUSIexon^I25237205I25237684^I.I+^I.Itranscript_id “g1.t1”; gene_id “g1”;$

sub(/\t+/,"\t")就是使用sub进行字符串替换。将多个tab替换成单个tab。
最后方案【！！！！】】】】】
$ awk ‘BEGIN {FS="[-:\t]";OFS="\t"}{if(/^#/){print $0}else{$6=$6+$2-1;$7=$7+$2-1;$2=$3="";sub(/\t+/,"\t");print}}’ target.gff > output.gff

bed文件坐标为一半开半闭区间[start, end),所以如果是[10,20)，实际上只提取了10,11，…19 这十个位点，对应ucsc上的即为染色体坐标的10-19位碱基。ucsc上染色体坐标也是从0开始。

例如有一fasta格式的文件

chr1

TCGAGA

对应bed文件的坐标应为

#chrome start end

chr1 0 5

用bedtools提取 CGAG 中间四个碱基，所需的bed输入文件应为[1,5)

#chrome start end

chr1 1 5

Bed文件和GFF文件最基本的信息就是染色体或Contig的ID或编号，然后就是DNA的正负链信息，接着就是在染色体上的起始和终止位置数值。
两种文件的区别在于，BED文件中起始坐标为0，结束坐标至少是1,； GFF中起始坐标是1而结束坐标至少是1。

从gff3中提取gene坐标等转为bed【没错 bed是需要gtf的start减一的】
awk ‘$3 == “gene”’ gencode.vM19.chr_patch_hapl_scaff.annotation.gff3 | awk ‘BEGIN{FS="\t|=|;";OFS="\t"}{print $1,$4,$5,$16}’ > genes.bed

去除非chr染色体的基因：

sed -n ‘/^chr/p’ genes.bed > genes.fine.bed

其他方法：

ref：RSeQC使用笔记 – 生信笔记

1.##gtf转化为bed：

cat reference/genome/hg19/gencode.v26lift37.annotation.gtf |awk ‘OFS="\t" {if($3==“transcript”) {print $1,$4-1,$5-1,$12,$6,$7}}’ |tr -d ‘";’ >hg19.bed

2.用ucsc的gtfToGenePred配合shell来将gtf转化为bed12：

gtfToGenePred -genePredExt -geneNameAsName2 ~/reference/genome/hg19/gencode.v26lift37.annotation.gtf gene.tmp
awk ‘{print $2"\t"$4"\t"$5"\t"$1"\t0\t"$3"\t"$6"\t"$7"\t0\t"$8"\t"$9"\t"$10}’ gene.tmp > hg19.bed12

这些列的数据，也可以用在circos上，进行基因组信息可视化。

求两个bed文件中描述的基因组区段的重叠区域
intersectBed可以用来求两个BED或者BAM文件中的overlap，overlap可以进行自定义是整个genome features的overlap还是局部。

bamToBed可以用来把比对产生的bam文件转化为Bed文件，可以配合coverageBed命令使用。

。。。。
1 基因组坐标转换常用软件
NCBI的Remap:支持BED、VCF、GFF、GTF等。
UCSC的LiftOver:支持BED,本地版本还支持GFF。
CrossMap:BAM, CRAM, SAM, BED, Wiggle, BigWig, GFF, GTF and VCF等。

2.1 安装
#安装(pip or conda)
conda install CrossMap
pip3 install CrossMap
chain文件作为转换坐标时候的输入文件，可以根据你的物种名在UCSC上下载(这里以小鼠mm9转换为mm10为例,下载mm9ToMm10.over.chain.gz文件即可)。
http://hgdownload.soe.ucsc.edu/downloads.html

#使用相对简单
CrossMap.py xx(什么文件格式:如bed) xx(chain文件) xx(原始坐标文件) xx(输出文件名)
#Note:由于我的bed文件莫名多出一列，这里先去掉一列，然后进行坐标转换。(注:bed文件里面第一行不能为列名，如果是请去掉)
awk ‘{$1="";print $0}’ gain_1K_sign_6.txt >gain_1K_sign_6.bed
CrossMap.py bed mm9ToMm10.over.chain.gz gain_1K_sign_6.bed gain_1K_sign_6_mm10.bed
#生成的结果map到的结果为xx.bed，同时也生成unmap的结果。

。。。。
生成特征文件！！chipseq 的overlap的基因！！！！

zcat gencode.v25.long_noncoding_RNAs.gtf.gz |perl -alne ‘{next unless $KaTeX parse error: Can't use function '\"' in math mode at position 28: …e" ;/gene_name \̲"̲(.*?)\";/; prin…$ F[0]\t $KaTeX parse error: Undefined control sequence: \t at position 5: F[3]\̲t̲$ F[4]\t$1" }’ >lncRNA.hg38.position
zcat gencode.v25.2wayconspseudos.gtf.gz |perl -alne ‘{next unless $KaTeX parse error: Can't use function '\"' in math mode at position 34: …t" ;/gene_name \̲"̲(.*?)\";/; prin…$ F[0]\t $KaTeX parse error: Undefined control sequence: \t at position 5: F[3]\̲t̲$ F[4]\t$1" }’ >pseudos.hg38.position
zcat gencode.v25.annotation.gtf.gz| grep protein_coding |perl -alne ‘{next unless $KaTeX parse error: Can't use function '\"' in math mode at position 28: …e" ;/gene_name \̲"̲(.*?)\";/; prin…$ F[0]\t $KaTeX parse error: Undefined control sequence: \t at position 5: F[3]\̲t̲$ F[4]\t$1" }’ >protein_coding.hg38.position

zcat gencode.v25.annotation.gtf.gz|perl -alne ‘{next unless $KaTeX parse error: Can't use function '\"' in math mode at position 28: …e" ;/gene_name \̲"̲(.*?)\";/; prin…$ F[0]\t $KaTeX parse error: Undefined control sequence: \t at position 5: F[3]\̲t̲$ F[4]\t$1" }’ >allGene.hg38.position

zcat gencode.v25lift37.annotation.gtf.gz | grep protein_coding |perl -alne ‘{next unless $KaTeX parse error: Can't use function '\"' in math mode at position 28: …e" ;/gene_name \̲"̲(.*?)\";/; prin…$ F[0]\t $KaTeX parse error: Undefined control sequence: \t at position 5: F[3]\̲t̲$ F[4]\t$1" }’ >protein_coding.hg19.position
zcat gencode.v25lift37.annotation.gtf.gz | perl -alne ‘{next unless $KaTeX parse error: Can't use function '\"' in math mode at position 28: …e" ;/gene_name \̲"̲(.*?)\";/; prin…$ F[0]\t $KaTeX parse error: Undefined control sequence: \t at position 5: F[3]\̲t̲$ F[4]\t$1" }’ >allGene.hg19.position

/data/zhangyong/yyp/yypold/ref/mm10/mm10gtf/allGene.mm10.position

cp /data/zhangyong/yyp/yypold/chips/lh724/sam/chenalign/chenpk/G*.xls /data/zhangyong/yyp/yypold/chips/lh724/sam/chenalign/overl/

cd /data/zhangyong/yyp/yypold/chips/lh724/sam/chenalign/overl
source activate rna
bedtools intersect -a G1-1_FKDL202583853_peaks.bed -b /data/zhangyong/yyp/yypold/ref/mm10/mm10gtf/allGene.mm10.position -wa -wb
| bedtools groupby -i - -g 1-3 -c 14 -o collapse >cs2【-c的第十行是-a文件的第14行，】

chr1 1000 2000 a1
chr1 5000 8000 a2

chr1 900 20002 b

bedtools intersect -a a.bed -b b.bed -wa -wb

cd /data/zhangyong/yyp/yypold/chips/lh724/sam/chenalign/overl
source activate rna
printf “chr1\t1000\t2000\ta1\nchr1\t5000\t80000\ta2” >a.bed
printf “chr1\t900\t20002\tb1” >b.bed
Error: Sorted input specified, but the file - has the following out of order record

printf “chr1\t900\t20002\ta3” >>a.bed

区间延长$ awk ‘BEGIN {FS="[-:\t]";OFS="\t"}{if(/^#/){print $0}else{$2=$2-2000;$3=$3+2000;print}}’ mG1.bed > 2kmG1.bed

bedtools sort -i 2kmG1.bed |bedtools merge > 2kmG1m.bed

bedtools intersect -a 2kmG1m.bed -b /data/zhangyong/yyp/yypold/ref/mm10/mm10gtf/allGene.mm10.position -wa -wb
| bedtools groupby -i - -g 1-3 -c 7 -o collapse >output

bedtools intersect -a 2kmG1m.bed -b /data/zhangyong/yyp/yypold/ref/mm10/mm10gtf/allGene.mm10.position -wa -wb|head

qq_39306047

关注

1
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
bed文件坐标的调整（awk的用法）

bed文件坐标的调整（awk的用法）。annotatePeak的输入也可以是GRanges对象，你如果用R做peak calling的话,直接就可以衔接上ChIPseeker了。require(ChIPseeker)f = getSampleFiles()[[4]]这里我们需要的是一个TxDb对象，这个TxDb就包含了我们需要的各种信息，ChIPseeker会把信息抽取出来，用于注释时使用。require(TxDb.Hsapiens.UCSC.hg19.knownGene)txdb
复制链接

扫一扫