how to retreive sequences given a start and end site.

So I have a gff that has information about MEF2 transcription factor binding sites. So given a start and end site, 19815641 - 19815654 on the + strand on ChrX, where exactly do I get the sequence from?

I have 1800 lines in the gff file, so I cant do it manually. I am looking for a R solution, so basically if something like the following function exists

getSequence(start, end, strand, chr)

The goal is to create a PWM so my next question is that once I've retrieved my sequences, how do I go about aligning them? what is the best software to align 1800 short sequences?

Edit: It seems like the

bedtools getfasta -fi reference.fasta -bed gff.file -fo output.fasta

is what I need, but whats the easiest way to download hg18 reference genome?

alignment sequence

• 2,471 views • link

updated 4.5 years ago by Ram

• written 11.7 years ago by Affan

http://hgdownload.cse.ucsc.edu/goldenPath/hg18/bigZips/

http://hgdownload.cse.ucsc.edu/goldenPath/hg18/chromosomes/

You can use samtools faidx to retrieve sequences from reference genome using coordinates.

faidx samtools faidx <ref.fasta> [region1 [...]]

Index reference sequence in the FASTA format or extract subsequence from indexed reference sequence. If no region is specified, faidx will index the file and create <ref.fasta>.fai on the disk. If regions are specified, the subsequences will be retrieved and printed to stdout in the FASTA format. The input file can be compressed in the RAZF format.

• 0 views • link

updated 4.5 years ago by Ram

• written 11.7 years ago by Ashutosh Pandey

0 answers

No answers yet.

More posts like this

Is there a way to remove certain chromosomes from a GRanges Object?
written 5.3 years ago by cthangav •

Hello, I am trying to do the tutorial here: https://raw.githack.com/bioFAM/MOFA2_tutorials/master/R_tutorials/10x_scRNA_scATAC.html But when I try to parse the seurat object using: motif.matrix <- CreateMotifMatrix( features = …
Bedtools maskFastaFromBed does not mask on the minus strand
written 6.5 years ago by Biostar •

Hello, I have an unmasked genome and its repeat file in bed/gff format. I tried to softmask the genome using maskFastaFromBed from bedtools (v2.26), but …
Retrieving sequences from Ensembl Archive
written 8.1 years ago by biostR •

Hi, I am looking for a way to retrieve DNA sequences from Ensembl May 2017 archive, based on coordinates. I thought using Biomart package would …
Bedtools with name and full header
written 10.4 years ago by bharata1803

Hello, So, I have a bed file and a fasta reference file. I want to generate a fasta file corresponding to the bed file. I …
bedtools getfasta does not give a proper fasta output
written 11.3 years ago by Rohit

I have a problem with the output from bedtools getfasta, it used to work before. My bed-file looks like this - C1_64090263 131 164 C1_64090263-1:33 …
calculating promoter region from TxDb for a negative strand gene
written 11.3 years ago by Affan

I am working on a R script that retrieves promoter regions for genes. My question isn't about the specific programming, but in a general biology …
How can I create a more accurate PWM?
written 11.6 years ago by Affan

I've created a PWM for the MEF2 transcription factor. I've done this using R. The MEF2 transcription factor binding sites were obtained from Riken database …
how to interpret scores from a PWM match
written 11.6 years ago by Affan

I created a PWM using true binding sites from Riken 4 database, hg 18. To test my PWM, I picked a random true binding site …
extracting parts of genome using R
written 11.6 years ago by Affan

Background: I downloaded TFBS data from the Riken database and created a PWM for the MEF2 transcription factor. Now I want to test the accuracy …
creating a PWM from TFBS data in a GFF file
written 11.7 years ago by Affan

So I downloaded the TFBS from the [Riken 4 database][1]. My task is to extract the binding sites for the transcription factor MEF2 and then …