Tag Info

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.

12

You can just set the fields you don't need to *: samtools view -h foo.bam | awk 'BEGIN{FS="\t"; OFS="\t"}{if($1~/^@/) {print$0} else {print "*", $2,$3, $4,$5, $6,$7, $8,$9, "*", "*"}}' | samtools view -bo smaller.bam - This will set each read's name to *, but you can still see where its mate maps with the PNEXT and RNEXT fields. The resulting BAM ...

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 ... 11 This can be done in R very easily from an indexed .bam file. Given single-end file for sample1. library(GenomicAlignments) library(rtracklayer) ## read in BAM file (use readGAlignmentPairs for paired-end files) gr <- readGAlignments('sample1.bam') ## convert to coverages gr.cov <- coverage(gr) ## export as bigWig export.bw(gr.cov,'sample1.bigwig') ... 11 If you have multi-line fasta files, as is very common, you can use these scripts1 to convert between fasta and tbl (sequence_name <TAB> sequence) format: FastaToTbl #!/usr/bin/awk -f { if (substr($1,1,1)==">") if (NR>1) printf "\n%s\t", substr($0,2,length($0)-1) else printf "%s\t", substr($0,2,... 10 There is a very simple BioPython solution, that is minimal, readable, and handles multi-line fasta: from Bio import SeqIO for record in SeqIO.parse('example.fa', 'fasta'): print('>{}\t{}'.format(record.description, record.seq)) 8 It's not a fasta file, but: > m sample1 sample2 sample3 sample4 fliI 1 1 1 1 patB_1 1 1 1 1 pgpA 1 1 1 1 osmB 1 1 1 1 cspA 1 0 1 1 > # Collapse to labeled strings > blah = apply(m, 2, function(x) paste(x, collapse=''... 8 assuming there is only one sequence line per record, use paste with two 'stdin' cat your.fasta | paste - - 7 You can use these commands: perl -pe 's/>(.*)/>\1\t/g; s/\n//g; s/>/\n>/g' sequences.fa | grep -v '^$' Explanation: Append a tab to every header line Join all lines Split the single obtained line by the '>' character Remove the empty line (the first line is empty due to the fact that '>' is the first character of the FASTA file)

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

NOTICE: I have altered my answer slightly from the original as I have turned the original script into a pip installable program (with tests) and have updated the links and code snippets accordingly. The essence of the answer is still exactly the same. This is something I have been meaning to get around to for a while, so thanks for the prompt. I have ...

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

I like bedtools getfasta. My typical option set is bedtools getfasta -fi <reference> -bed <gff_file> -name -s. Be aware of the -s to make sure you are pulling the correct strand. I like bedtools because it is a versatile tool overall for handling bed, gff and vcf file manipulations. # bedtools getfasta Tool: bedtools getfasta (aka ...

6

According to the documentation (?Biostrings::DNAStringSet): width(x): A vector of non-negative integers containing the number of letters for each element in x. Note that width(x) is also defined for a character vector with no NAs and is equivalent to nchar(x, type="bytes"). names(x): NULL or a character vector of ...

5

I'm not aware of any pre-made program to do this, so I wrote one for you. This will take a BAM file with any ordering and produce properly ordered gzipped fastq files with the filtering as you requested. Internally, this iterates over all of the entries in the BAM file (ignoring secondary/supplemental entries and those where both mates map to your filter ...

5

The following bit of python code should work: #!/usr/bin/env python import sys lastTranscript = [None, None, None, []] # ID, chrom, strand, [(start, end, score), ...] def getID(s): """Parse out the ID attribute""" s = s.split(";") for k in s: if k.startswith("ID="): return k[3:] return None def dumpLastTranscript(): ...

5

You can do this in Hail. Here's the rough code to do it (0.1 versions). Setup: from hail import * hc = HailContext() Import the .gen file. VCF works too: dataset = hc.import_gen( 'src/test/resources/example.gen', 'src/test/resources/example.sample') Remap the genotype schema and export to VCF: dataset.annotate_genotypes_expr('g = {GT: g.call()...

5

A very useful tool for this kind of data manipulation is bioawk: $bioawk -c fastx '{print ">"$name" "$comment"\t"$seq}' test.fa >sample 1 gene 1 atgc >sample 1 gene 2 atgc >sample 2 gene 1 atgc bioawk is based on awk, with added parsing capabilities. Here, we tell that the format is fasta or fastq with -c fastx, and this makes the $... 5 Finally, I found an alternative to the SRA translation: a link that works! For those of you interested in knowing how to download FastA files from NCBI using an accession number, try the following link: https://www.ncbi.nlm.nih.gov/search/api/sequence/${accession}/?report=fasta Using wget to download the accession used as example: wget -O NC_001416.1....

5

A gff is a file of annotation. It generally doesn't include sequence information, so you can not inter-convert. But the file on your page has a fasta entry at the bottom.. I'd just grab the bottom half of it, and if the sequence length looks right, that's likely what you want.

5

You can use GFF utilities to get a fasta sequence from a GFF file using a reference genome file.

5

Atom reordering is a common problem in compchemistry. Rdkit (a python package) can do this, but it is limited by the formats it can read and mol2 files are a bit hit or miss. It works really well with SMILES, SMARTS and mol (sdf) files. But the writing may cause problems with Brookhaven pdb and mol2 files. So the formats you have are both problematic. One ...

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

There's no equivalent to the wiggle header in bigWig (or bigBed) files, which is why UCSC uses the file name. This is actually the reason for the track line stuff that you linked to, since you can then specify a name and just point to where the bigWig (or other format) file is on the internet. BTW, you can certainly convert your bigWig to wiggle, add the ...

4

Apparently convert_trace does not do a good parameter checking and silently sets ZTR as default if it does not recognise a valid output format. I don't know of any freely available command-line converters in that specific direction. Even at the height of Sanger sequencing (late 90s, early 2000s), pipelines usually tried to reduce complexity and footprint by ...

4

Just for reference, there is a not very in depth, but first-hand reasoning (Jim Kent) given here about why bedGraphToBigWig does not support streaming http://genome.soe.ucsc.narkive.com/2S1Z3VpG/bedgraphtobigwig-reading-in-from-stdin Edit: as noted wigToBigWig allows bedgraph streaming input via stdin but takes (possibly much) more memory than ...

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