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/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+