extract transcript boundaries

December 14, 2018 · View on GitHub

requirement bed tools

BIN='/home/hirak/bedtools2/bin'

Gencode

gencode.v29.chr_patch_hapl_scaff.annotation.gtf

GTF_FILE="gencode.v29.chr_patch_hapl_scaff.annotation.gtf"

extract transcript boundaries

cat GTF_FILE | awk 'BEGIN{OFS="\t";} \3=="transcript" {print $1,$4-1,$5,$12}' | tr -d """ | tr -d ";" | $BIN/sortBed > gencode_transcript_intervals.bed

merge exon boundaris

cat GTF_FILE | awk 'BEGIN{OFS="\t";} \3=="exon" {print $1,$4-1,$5,$12}' | tr -d """ | tr -d ";" | BIN/sortBedBIN/sortBed | BIN/mergeBed -i - -c 4 -o collapse > gencode_exon_merged.bed

extract introns per transcript

$BIN/subtractBed -a gencode_transcript_intervals.bed -b gencode_exon_merged.bed -nonamecheck > intron_transcript.bed

Group by by awk

#awk -F "\t" '{a[$1"\t"$2"\t"$3]=a[$1"\t"$2"\t"$3]?a[$1"\t"$2"\t"$3] OFS $4:$4 } END {for (i in a) print i FS a[i]}' OFS=, intron_transcript.bed | /mnt/scratch1/hirak/bedtools2/bin/sortBed -i -

Or by bedtools groupby

sort -k1,1 -k2,2 -k3,3 intron_transcript.bed | $BIN/groupBy -i - -g 1,2,3 -c 4 -o collapse > hg_intron_transcript_collapsed.bed

make fasta by using this bed

$BIN/bedtools getfasta -fi hg_genome.fasta -bed intron_transcript.bed -fo hg_intron_transcript.fasta -name+