Hot answers tagged

23

You're the second person I have ever seen using NCBI "chromosome names" (they're more like supercontig IDs). Normally I would point you to a resource providing mappings between chromosome names, but since no one has added NCBI names (yet, maybe I'll add them now) you're currently out of luck there. Anyway, the quickest way to do what you want is to samtools ...


14

You can do this easily with bioawk, which is a version of awk with added features facilitating bioinformatics: bioawk -c fastx '{print $name"\t0\t"length($seq)}' test.fa -c fastx tells the program that the data should be parsed as fasta or fastq format. This makes the $name and $seq variables available in the awk commands.


13

Via Gencode and BEDOPS convert2bed: $ wget -qO- ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gff3.gz \ | gunzip --stdout - \ | awk '$3 == "gene"' - \ | convert2bed -i gff - \ > genes.bed You can modify the awk statement to get exons, by replacing gene with exon. BEDOPS: https://github.com/...


11

It's good practice to have your FASTA indexed, so you can leverage the .fai you are likely to already have. If not, you can just generate the index with samtools and use some awk to make your BED: samtools faidx $fasta awk 'BEGIN {FS="\t"}; {print $1 FS "0" FS $2}' $fasta.fai > $fasta.bed This will maintain tab separation but you can drop the BEGIN ...


8

I am unaware of any "official" or gold-standard UTR annotations in S. cerevisiae. One option is to use the annotations from the TIF-Seq publication (Pelechano et al. 2013). The GSE39128_tsedall.txt.gz file contains the major isoforms identified. It would be up to you to computational associate each transcript with a given gene. It is also up to you to ...


7

You can store it in data/ subfolder. If one of your functions need this bed file, you can include importing it in your function. Here are some examples how to get the data from data/ subfolder. There are probably more ways to do this, this is just one way.


6

Galaxy has API and API-consuming libraries (such as BioBlend) that will allow you to interactively script against it without opening the graphical interface at all. However you can also take almost any tool out of Galaxy and use it independently since everything is open source. The converter you mentioned is available as a Python script here and the tool '...


6

You could adapt this awk one-liner. Note that it assumes that sequence IDs are not longer than 100 characters and that there is no description following the sequence ID on the header line. cat myseqs.fasta | awk '$0 ~ ">" {print c; c=0;printf substr($0,2,100) "\t0\t"; } $0 !~ ">" {c+=length($0);} END { print c; }' Otherwise, any Bio* library (Perl, ...


6

According to the SAM specification, the 3rd field of a SAM line (RNAME) is: RNAME: Reference sequence NAME of the alignment. If @SQ header lines are present, RNAME (if not ‘*’) must be present in one of the SQ-SN tag. An unmapped segment without coordinate has a ‘*’ at this field. However, an unmapped segment may also have an ordinary coordinate ...


6

As benn said you can store them in data/ subfolder of your package. There is no size limit (AFAIK) but if you plan to publish your package in some repositories they might have limits or other restrictions. Bioconductor requires to divide the data required for a software package at certain threshold, but for small amount of data it is fine. You can find ...


6

I found the following two files in https://downloads.yeastgenome.org/sequence/S288C_reference/: SGD_all_ORFs_3prime_UTRs.fsa SGD_all_ORFs_5prime_UTRs.fsa According to the README files in the same directory, these are (the README for the 5' file is equivalent): SGD_all_ORFs_3prime_UTRs.README Information about the SGD_all_ORFs_3prime_UTRs.fsa file. ...


6

awk 'BEGIN{OFS="\t"}{if($1==lchrom && $4==lname && $2 == lend) {lend = $3}else{if(lchrom) {print lchrom, lstart, lend, lname;}; lchrom=$1; lstart=$2; lend=$3; lname=$4}}END{print lchrom, lstart, lend, lname}' foo.bed > new.bed That produces: chr1 0 10000 LOC101929784::chr1:13201-13800::E1 chr1 10000 10600 LOC101929784::chr1:...


5

You're reinventing bedtools intersect (or bedops), for which there's already a convenient python module: from pybedtools import BedTool s3 = BedTool('s3.bed') s4 = BedTool('s4.bed') print(s4.intersect(s3, wa=True, wb=True, F=1)) The wb=True is equivalent to -wb with bedtools intersect on the command line. Similarly, F=1 is the same as -F 1.


5

I've written a handful of programs from scratch to simulate mutations and variations in real or simulated sequences. The trick has always been to sort the variants by genomic coordinate, apply the variant with the largest coordinate first, then apply the variant with the second largest coordinate, all the way down to the variant with the smallest coordinate....


5

Although you don't mention it, I'm guessing you're using bedtools v2.26.0. Version 2.26.0 of groupBy has a bug in it, which you've encountered (it was fixed shortly after release, so you'll either have to use a version before the bug was introduced, or compile the current source code yourself from https://github.com/arq5x/bedtools2) v2.26.0: local10:~/...


5

There are multiple ways to go about it. On the command line you can make a 1 line BED file: chr1 11868 12227 And then bedtools intersect with it. In R, you could load your original BED file and use GenomicRanges: library("GenomicRanges") bed = read.delim("foo.bed", header=F) # Rename this # N.B., BED files use 0-based coordinates, I've switched to 1-...


5

# if you have seqtk installed, skip the following two lines git clone https://github.com/lh3/seqtk cd seqtk && make # the main command line ./seqtk cutN -gp10000000 -n1 hg38.fa > hg38-N.bed Option -n sets the min stretch length. Just use -p as is. It is a bit complicated to explain what it is doing.


5

Here's a way to use BEDOPS, which was designed to work fast by using sorted input. Other tools now use sorting to accomplish similar performance benefits. Convert GTF annotations to a sorted BED file of genes: $ awk '($3=="gene")' annotations.gtf | gtf2bed - > genes.bed Sort your intervals, if unsorted: $ sort-bed intervals.unsorted.bed > intervals....


5

In one line, using bedtools zcat Homo_sapiens.GRCh38.93.gtf.gz \ | awk '$3=="gene"' \ | bedtools slop -b 10000 -g contigs.tsv -i - \ | bedtools intersect -u -a intervals.bed -b - This first takes the genes and filters them for the third column being gene. This is important because all GTF lines contain the word gene as gene_id is a manditory attribute ...


4

You could do this with the CGAT toolkit: cgat bed2bed --method=merge --merge-by-name -I bed_with_gene_ids.bed Installing such a massive package might be overkill for this task though.


4

To answer the question as asked, for people googling. For BED6, in python: #contigs.tsv contians chromosome names and lengths in two columns for line in open("contigs.tsv"): fields = line.strip().split("\t") print fields[0], ".", "contig","1",str(fields[1]), ".", "+", ".", "ID=%s" % fields[0] for line in open("my_bed_file.bed"): fields = line....


4

We have many excellent answers! This will be an excellent reference for future users. I found what exactly what I was asking in my question: https://www.biostars.org/p/191052/ $ pip install pyfaidx $ faidx --transform bed test.fasta > test.bed This is the one-line command I was asking. The other answers also work, but I want to accept my own ...


4

The "right" solution would be realignment, but that's expensive and most of us would not go that route. My preferred solution would be to convert the bed file, as opposed to the bam. Here's why: 1) Reheadering the bam means that you may have reads aligned to contigs without a corresponding entry in UCSC (see Devon's list for the mappings). This is a problem ...


4

Here is an approach with BioPython. The with statement ensures both the input and output file handles are closed and a lazy approach is taken so that only a single fasta record is held in memory at a time, rather than reading the whole file into memory, which is a bad idea for large input files. The solution makes no assumptions about the sequence ID lengths ...


4

Inspired by this answer to a related question on read length distributions, you could do this with Biopython: from Bio.SeqIO import parse with open("regions.bed", "w") as bed: for record in parse("regions.fasta", "fasta"): print(record.id, 0, len(record.seq), sep="\t", file=bed)


4

If you look for a program which would randomly introduce SNPs + short indels and then would save everything into a VCF file, DWGsim or Mason Variator could be a good choice. Then you can create a corresponding Chain file using bcftools consensus -c and transform various formats between these two coordinate systems using CrossMap.


4

I quite like BED and GFF3 (I don't like GTF/GFF2, though). As text-based formats, I don't think they leave us much room for improvement. Anyway, if you want a new format, here is one. The following is a hybrid between GFF3 and BED. It is a TAB-delimited text-based format with the following fields: chr, required start (0-based), required end, required strand ...


4

The problem is that GFF, fundamentally, is a relational format: it provides tags that relate individual rows via one-to-many relationships (e.g. gene–exon). This indirectly highlights the second complication: individual rows have different types, and therefore store different attributes in the 9th column. Over the last few decades (!), we have accumulated a ...


4

Presuming we consider "human readable", "easily parsable", and "quickly queryable" to be objectively good qualities (and if not, I worry for the future): Text-based: It's absurdly common to want to use grep or awk on annotations. Sure, one could make variants of these that are binary-format aware, but why reinvent the wheel. Of course text files don't ...


4

Figured it out! For anyone who finds this thread in the future and is wondering what is going on: HOMER, when using the -size given option, normalizes the y-axis to the number of basepairs in the intervals -- in this case, the length of the gene body. To get the "raw" values, as I wanted, you need to multiply by the length of the gene body to "undo" the ...


Only top voted, non community-wiki answers of a minimum length are eligible