How can I get the sequences of the first 20nt of introns and the last 250nt of introns of hg38 genome?

Question

[xiaoleiusc]

Dear Biostars users,

I have the hg38 genome fasta file and a bed file of all the introns (GENCODE V44) of the hg38 genome. I would like to get two fasta files: 1) a fasta file of the first 20nt of all the introns; and 2) a fasta file of the last 250nt of all the introns.

What tools (and command lines) should I use to generate these two output fasta files?

Thanks ahead,

Xiao

Answer 1

[harold.smith.tarheel]

You can use Bedtools flank to generate the BED file of the desired intron intervals, then getfasta to output the sequences.

[xiaoleiusc]

The example code is bedtools flank [OPTIONS] -i <BED/GFF/VCF> -g <GENOME> [-b or (-l and -r)], in my case, can I use -l -20 and -r -250 ? I wonder if I can use negative values (e.g. -20 or -250 in my case) followed the -l or -r. Thanks.

[GenoMax]

-l  The number of base pairs to subtract from the start coordinate. Integer.

You are covered on the left end. May need to create a BED with one base intervals with start and end of the introns. One entry for start and one for end.

Answer 2

[Juke34]

Another solution is to use AGAT Either you use the gff/gtf annotation as input or convert your bed to gff with agat_convert_bed2gff. Then you can use agat_sp_extract_sequences using parameter to extract flank regions : https://agat.readthedocs.io/en/latest/tools/agat_sp_extract_sequences.html