cannin/biostars_qa huggingface dataset:
Question: Hi, I want to convert my text file to a GTF, GFF3 or BED file. My file is organized like that.
This is a list of operons downloaded from RegulonDB, but I need to convert it in one of this format in order to visualize them in IGV.

Column Names:
1) OPERON_ID
2) OPERON_NAME
3) FIRSTGENEPOSLEFT
4) LASTGENEPOSRIGHT
5) REGULATIONPOSLEFT
6) REGULATIONPOSRIGHT
7) OPERON_STRAND
8) KEY_ID_ORG

ECK120007527	ycaM	946452	947882	946159	947882	forward	ECK12
ECK120014360	astCADBE	1823979	1830006	1823979	1830351	reverse	ECK12
ECK120014361	pyrH	191855	192580	191718	192580	forward	ECK12
ECK120014362	nrdHIEF	2798745	2802483	2798457	2802483	forward	ECK12
ECK120014363	cpxP	4103843	4104343	4103749	4104402	forward	ECK12

I do not know if there is a way to directly convert my file in one of these formats. Otherwise , what is the command to get the correct columns ?
Thank you.

Answer: Hi bonardi.franck,

All the files you mentioned are text-files with a certain formal definition. See here the [bed][1], the [gff3][2], and [gtf][3] descriptions.

For your purpose, it seems to me it would be best to parse your data into a bed format. You can use FIRSTGENEPOSLEFT and LASTGENEPOSRIGHT to define your region (fields 2 and 3 ) and the REGULATIONPOSLEFT and REGULATIONPOSRIGHT as thick start, resp end. For the name, you can choose the OPERON_ID, the OPERON_NAME or a mixture of both. The OPERON_STRAND needs to be converted into + or -. Still you need the chromosome name (I assume KEY_ID_ORG is this here).

You can use for instance awk or perl to reformat each line of your textfile, that you have the necessary structure. It looks something like `awk '{if(match($7,/forward/)){s="+"}else{s="-"} ; print $8"\t"$3"\t"$4"\t"$1"\t"s"\t"$5"\t"$6 }' mytext.txt > myoperons.bed`

Cheers,

Michael

[1]: http://bedtools.readthedocs.org/en/latest/content/general-usage.html
[2]: http://www.sequenceontology.org/gff3.shtml
[3]: http://mblab.wustl.edu/GTF2.html
Question: Hey guys,

I am going to aligning reads from GSE35641 (http://www.ncbi.nlm.nih.gov/Traces/sra/?view=search_seq_name&exp=SRX119208&run=&m=search&s=seq), on ecoli reference genome (NC_000913). There are 6 sra files in GEO which the least size is more than 2G. Do you know any way to download the files faster? I need many hours to download a file of them in fasta format

Thank you biostars

Answer: Look at using [aspera from the command line][1] or, if you are using R, as described in sections 3.5 and 3.6 of [the SRAdb vignette][2].

[1]: http://www.ncbi.nlm.nih.gov/books/NBK158899/
[2]: http://www.bioconductor.org/packages/release/bioc/vignettes/SRAdb/inst/doc/SRAdb.pdf
Question: Hi,

I have two files with about 4 fasta sequences each in the PHYLIP format for Da/Ds (Ka/Ks) comparison. After a search on the net, I found that the best software to do so is PAML and it requires a control file and 2 input files. For the tree input, I used PHYML which accepts PHYLIP file inputs and gives NEWICK format outputs. However I am unsure as to how to get the sequence file. I tried the sequential format but it gives me the "Error in sequence data file" error.

[This][1] is my control file. [Here][2] is my sequence file and [here][3] is my tree file.

Could I get some advice on what is wrong with my input file please. Also, this is just one set of the sequences which I wanted to compare. How do I compare two sets of sequences for the Da/Ds values. I know that this might be basic so thank you for your patience.

[1]: https://www.dropbox.com/s/dupmdmufapmfzdk/codeml.ctl?dl=0
[2]: https://www.dropbox.com/s/mwr05pxvlzvk9ea/tm-lrr.txt?dl=0
[3]: https://www.dropbox.com/s/e6hi4pirbr8goho/tmphy_phy_phyml_tree.txt?dl=0

Answer: I don't think your sequence file is in correct format and does not look like a multiple alignment. You need a multiple sequence alignment as input. You can use the online tools (either Readseq or Built-in converter) at [Phylogeny.fr][1] website to check if your sequence file is in correct format (by trying to convert to other format).

The sequences in your file are also not of same length. If the 'cleandata' option is set to 1 in the .ctl file (which is the default), PAML will remove any sites with gaps in at least one sequence.

From the [PAML User Guide][2] (Page 11):

The way that ambiguity characters and alignment gaps are treated in baseml and codeml depends on the variable `cleandata` in the control file. In the maximum likelihood analysis, sites at which at least one sequence involves an ambiguity character are removed from all sequences before analysis if `cleandata` = 1, while if `cleandata` = 0, both ambiguity characters and alignment gaps are treated as ambiguity characters. In the pairwise distance calculation (the lower-diagonal distance matrix in the output), `cleandata` = 1 means "complete deletion", with all sites involving ambiguity characters and alignment gaps removed from all sequences, while `cleandata` = 0 means "pairwise deletion", with only sites which have missing characters in the pair removed.

[1]: http://phylogeny.lirmm.fr/phylo_cgi/programs.cgi
[2]: http://abacus.gene.ucl.ac.uk/software/pamlDOC.pdf
Question: The lastz command will look like:

lastz reference.fasta sample.fasta --ambiguous=IUPAC --format=GENERAL > output.fasta

I need to run lastz once for every sample file, because each sample file has to be compared to its own individual reference (I have one reference for each sample).

I'm using find to get all of my 'sample' files, and will do the same for the reference files.

find . -type f -name 'sample.fasta' -not -path "*/temp/*"

The question is, how do I associate each sample file with each reference file, taking as input the output of the find command? Or is there a better way?

Answer: $ tree .
.
├── a
│ ├── reference.fa
│ └── sample.fa
├── b
│ ├── reference.fa
│ └── sample.fa
└── c
├── reference.fa
└── sample.fa

3 directories, 6 files

using parallel (remove dry-run option to execute the code):

$ ls -d */ | parallel --dry-run 'lastz ./{}reference.fa ./{}sample.fa --ambiguous=IUPAC --format=GENERAL > ./{}output.fa'

commands to be executed (i.e dry-run of code):

lastz ./a/reference.fa ./a/sample.fa --ambiguous=IUPAC --format=GENERAL > ./a/output.fa
lastz ./b/reference.fa ./b/sample.fa --ambiguous=IUPAC --format=GENERAL > ./b/output.fa
lastz ./c/reference.fa ./c/sample.fa --ambiguous=IUPAC --format=GENERAL > ./c/output.fa
Question: Hi,
I am generating bcf files from vcf. When indexing using bcftools the default indiexing is .csi however I also have tabix (generates .tbi index files) installed on my system hence I can opt for either. I'd just like to have a brief comparison between the two indexing formats in order to know which one to choose.

Thank you,

Answer: Hello ThePlaintiff!

Questions similar to yours can already be found at:

samtools 1.0 'index' : CSI index vs BAI index ?

We have closed your question to allow us to keep similar content in the same thread.

If you disagree with this please tell us why in a reply below. We'll be happy to talk about it.

Cheers!
Question: I have 21 bacterial isolates from termite guts and I want to determine the completness of the genomes. Which tool between BUSCO and CheckM is best to do this?

Answer: I use CheckM and haven't tried BUSCO in couple of years, so take that into account when considering my advice.

My impression was that BUSCO was too stringent in how it defined completeness and that they didn't use enough reference groups. That might have changed with version 4 of BUSCO based on the increased number of taxonomic groups, but I have no first-hand experience with the latest release. I think CheckM is easy to install and its results are straightforward, though it does require large memory. It is also my impression that CheckM is more or less the standard for prokaryotes (does not work with eukaryotes), while BUSCO may be more useful for eukaryotes even though it works with both. In the long run, BUSCO might overtake CheckM because it is broader in scope and appears to be actively maintained.
Question: I am in the typical situation that I need a resequencing pipeline (i.e., FastQC, read preprocessing, FastQC again, alignment with BWA, variant calling). I need to fulfill both the requirements of having a stable pipeline with stable tools for the standard stuff (e.g., both "single-donor WES variant calling", "trio WES variant calling", but also "tumor/normal WES variant calling with somatic filtration") but I sometimes need more specialized functionality or more extensive downstream analysis.

I want to use Docker for isolating my tools against the uncertain, changing, and sadly oftentimes unversioned world of Bioinformatics software (I'm looking at you, vt and vcflib, but I'm still very grateful that you are around). What would be your recommendation for a best practice here:

- one Docker image for everything, adding tools as I go
- one Docker image for each pipeline step (e.g. combining BWA-MEM, samtools, samblaster for the alignment so I can use piping in a front-end script)
- one Docker image for the standard stuff, then maybe some images for each additional step.

Does anyone know of a person/organization that has published their Dockerized pipeline stuff in a Blog post or elsewhere that goes beyond toy examples or "here is a Dockerfile for the tool that I wrote/published"?

Cheers!

Answer: Illumina Base Space actually uses it to run their cloud based pipeline with Docker tools. http://basespace.illumina.com They keep it down to tools level.

Docker was initially built and optimized to run "one" process efficiently. Therefore it is advisable to keep it down to tools level. Also, if you think about pipeline building, you probably want to have a Lego-like tool box and combine different workflows with different - often the same - tools.

As Illumina solved this, they have in each Docker container an input folder and an output folder, in Basespace you can actually publish your own tools as Docker container.
Question: Hi, I have the gtf file for mouse, and I wanted to create bed files for the 5' and 3' UTRs using this gtf file ( I know there is method using the UCSC Table website, but I want to generate using the GTF file itself).

As far as I know, information about the UTR regions is not explicitly stated in the gtf file. So, what would be a good strategy to extract these bed files from the gtf file?

Thanks!

Answer: Another [option][1].

python extract_transcript_regions.py -i genes.gtf -o genes --gtf

now convert this blockbed (bed12) to bed6

cat genes_3utr.bed | bed12ToBed6 -i stdin -n > genes_3utr_bed6.bed
cat genes_5utr.bed | bed12ToBed6 -i stdin -n > genes_5utr_bed6.bed

[1]: https://github.com/stephenfloor/extract-transcript-regions
Question: Hello

I want to parallelize a fastQC computation of a large fastq file.

I know that fastQC can be run on several files in parallel using the -t option, but I'm not sure if this will yield a correct result if I will split the file into several smaller files.

Is there any way to this?

Thanks in advance
Dolev Rahat

Answer: If you split the files into small parts then FastQC will run on each separately and the results are not combined.

As to an answer to your question - the problem (as always) is more complicated than it looks. In a perfect world splitting the file into random subsets would produce the similar results. But splitting into random subsets is not all that easy. In addition the way FastQC works is that it some operations will only collect information from the first 200K or so reads but will then track that information throughout the file.

Often and if the data is not skewed in particular ways you can get by analyzing a small subset of the file.

In other cases (as I have painfully learned myself) subsetting a dataset can produce reports that do not even remotely resemble the original data. But that had nothing to do with the way FastQC works.

Long story short I would run FastQC on all data then a subset of it - check if the library prep and the subject under study supports subsetting.
Question: I have a combined fasta file with all my sequences. I want to print the ID lines and their DNA sequence that end in `.4`

So far I have `awk /'>*.4'/ {getline;print} combined.fasta` Which prints the sequences that I want, how do I get the ID lines as well?

Answer: Since you requested bash only:

#!/bin/bash

string="$2"

while read line ; do
if [[ ${line:0:1} == ">" ]] ; then
header="$line"
else
seq="$line"
if [[ "$header" == *"$string" ]]; then
echo -e "$header""\n""$seq"
fi
fi
done < $1

Put it in a script file and run it like:

$ bash parseheaders.sh seqs.fasta .4

Some of my own test data:

$ cat seqs.fasta
>tpg|Magnaporthiopsis_incrustans|JF414846
ACTGTAGTAGCTACGATCGATCAGATGATCACGTAGCATCGATCGATCATCGACTAGTAGATCACTCGACATAGATCCACATCAATAGATCATCATCATCATAATCGATCACTAGCAGC
>tpg|Pyricularia_pennisetigena|AB818016
GCAAGNTTCATGACGATGTAGAATGGCTTATCGAAGGGAGCAGGCCAGGGATTGAGGTCCGTCTCACGGGTTGGCTTCACTCCCCCACTGCCAGCCCTCTTGCTGCAACTCCACCAGAA
>tpg|Inocybe_sororia|EU525947
AACCANGCCGCGACGGCGGTGCGATCGGGAAACGCGGCGGTGGCGGAGGAATCGGCCATCCTTCACCATATCGGCCAAGGATTGTGGTTCCTGTAGGGCTCGCGCAGCCCAGGACGCGC

$ bash parseheaders.sh seqs.fasta 947
>tpg|Inocybe_sororia|EU525947
AACCANGCCGCGACGGCGGTGCGATCGGGAAACGCGGCGGTGGCGGAGGAATCGGCCATCCTTCACCATATCGGCCAAGGATTGTGGTTCCTGTAGGGCTCGCGCAGCCCAGGACGCGC
Question: Hi group,

Just a quick question, how to I properly use the fastq command in samtools 1.3? According to the documentation I should be doing this:

samtools fastq -On -s ion36_nohost.fastq ion36_nohost.bam

but the result is this:

[M::bam2fq_mainloop_singletontrack] discarded 0 singletons

[M::bam2fq_mainloop_singletontrack] processed 0 reads

I'm trying to get the unmapped reads. My bam file does contain the reads in question but I'm doing something wrong since they switched to `fastq` from `bam2fq`

Any idea?

Answer: I must have misunderstood the docs becasue this was at the top and works

samtools fastq input.bam > output.fastq

but the detail doc is different here:

> samtools fastq [options] in.bam samtools fasta [options] in.bam
>
> Converts a BAM or CRAM into either FASTQ or FASTA format depending on
> the command invoked.
>
> OPTIONS:
>
> -n By default, either '/1' or '/2' is added to the end of read names where the corresponding BAM_READ1 or BAM_READ2 flag is set. Using -n
> causes read names to be left as they are.
>
> -O Use quality values from OQ tags in preference to standard quality string if available.
>
> -s FILE Write singleton reads in FASTQ format to FILE instead of outputting them.
>
> -t Copy RG, BC and QT tags to the FASTQ header line, if they exist.
>
> -1 FILE Write reads with the BAM_READ1 flag set to FILE instead of outputting them.
>
> -2 FILE Write reads with the BAM_READ2 flag set to FILE instead of outputting them.
>
> -0 FILE Write reads with both or neither of the BAM_READ1 and BAM_READ2 flags set to FILE instead of outputting them.
>
> -f INT Only output alignments with all bits set in INT present in the FLAG field. INT can be specified in hex by beginning with `0x' (i.e.
> /^0x[0-9A-F]+/) or in octal by beginning with `0' (i.e. /^0[0-7]+/)
> [0].
>
> -F INT Do not output alignments with any bits set in INT present in the FLAG field. INT can be specified in hex by beginning with `0x'
> (i.e. /^0x[0-9A-F]+/) or in octal by beginning with `0' (i.e.
> /^0[0-7]+/) [0].
Question: HI,

I need some help to start working with SAMtools. As a beginner with biostatistics background, it is my first time that I'm using Linux Ubuntu (I have installed it on my windows by making a machine on VMware for illustration purpose). Now I'm going to know how to work with SAMtools to carry out basic commands such as:

1- importing bam/sam files into samtools
2- Viewing the data
3- Sorting bam files
4. Calling variants

I have already searched on the internet and have found some tutorials such as http://www.htslib.org/doc/samtools.html , but I don't have any idea what those commands are!

could anybody please tell me where I can find a step by step tutorial for a beginner who even doesn't know how to read such commands?

If I am right, I should run samtools on terminal... by using the commands : "sudo apt-get install samtools" and then writing "samtools", I think I have run samtools but how should I write the commands in samtools on terminal?

I would also be grateful if someone tells me where I can find some low sized bam files to test samtools.

Thanks a lot!

Answer: Spend some time here learning the basics of [UNIX][1].

You can also take a look at the online materials for a NGS data analysis course that can be found [here][2]. This would give you an overall idea of steps you would need to learn and the materials for that.

[1]: http://korflab.ucdavis.edu/Unix_and_Perl/current.html#part1
[2]: http://angus.readthedocs.org/en/2014/
Question: Dear All,

I am searching for human 28S, 18S, 5.8S and 5S rRNA (RNA28S, RNA18S, RNA5-8S and RNA5S genes) and also 12S and 16S mitochondrial rRNA (MT-RNR1 and MT-RNR2) in human_genes.gtf file from UCSC.

I am unable to find these genes(RNA28S, RNA18S, RNA5-8S, RNA5S, MT-RNR1 and MT-RNR2 genes) in `human_genes.gtf` file. I am looking at `gene_id` or `gene_name` in gtf file.

How do I check these genes in gtf file or in any other resources.

Is it saved in some other name?

Answer: You can find 5S, 5.8S and MT_rRNA in the Ensembl/Gencode files. The others aren't in the reference genome (they're on an unplaced contig). If you need the sequences then you can get them from NCBI (the ones not in the reference genome come from 45S).
Question: Hi,

I want to extract reads from bam file using Samtools view from multiple regions.The regions I have are present in the VCF file,which I want to use.So there are not a small number of regions.Is there any way to do this in Samtools view.

for e.g inputting a file in some format in Samtools view?

Thanks

Answer: Hello Ron!

Questions similar to yours can already be found at:

Samtools And Region List

We have closed your question to allow us to keep similar content in the same thread.

If you disagree with this please tell us why in a reply below. We'll be happy to talk about it.

Cheers!

PS: duplicate of "Samtools And Region List"
Question: Hi all ...

I'm using Biopython 3.4.1 with Blast 2.2.29 and I'm trying to get the query coverage of some blast search beside the default report info in tabular format using `-outfmt 6` as following :

comline=NcbiblastxCommandline(cmd='blastn',query="A.fasta",db="B.db",evalue=0.001,outfmt='6 qcovs',out="result.txt")

But it throws the following error:

Error: Too many positional arguments (1), the offending value: qcovs

I know it seems a very basic question ,but it's seriously driving me crazy ,since I spent two days trying to figure out the solution for this error and nothing seems to fix it !!

Some threads have addressed this issue and suggested to rewrite the command by hand instead of copying and pasting ... I did rewrite it by hand like thousand times and the error pertains .

Also, I have tried to run it from the Command Prompt cmd and again it threw the same error!!

Is there anyway to get around this problem?

Many thanks in advance

Answer: Did you try `outfmt='"6 qcovs"'`? The command string needs to literally contain these quotes, such that the command line looks like:

blastn -outfmt "6 qcovs" ...
Question: Hi Every one

I am a bit stuck at this error from Michigan Imputation Server, does any one know what this error means and any hint to fix?

> Error: More than 100 obvious strand flips have been detected. Please
> check strand. Imputation cannot be started!

Statistics:
Alternative allele frequency > 0.5 sites: 0
Reference Overlap: 99.80%
Match: 426,485
Allele switch: 66,772
Strand flip: 23,263
Strand flip and allele switch: 13,805
A/T, C/G genotypes: 642
Filtered sites:
Filter flag set: 0
Invalid alleles: 0
Duplicated sites: 286
NonSNP sites: 0
Monomorphic sites: 0
Allele mismatch: 465
SNPs call rate < 90%: 0

Thanks

Answer: beagle provided a java tool called `conform-gt` https://faculty.washington.edu/browning/conform-gt.html

a program for making alleles in a VCF file to be consistent with a reference VCF file, see an example:

```
java -jar conform-gt.24May16.cee.jar gt=input.vcf.gz chrom=$i ref=chr22.1kg.phase3.v5a.EUR.vcf.gz out=output.chr22.beagle.vcf.gz
```
Question: I am attempting to re-create this figure from this link http://figshare.com/articles/_Enrichment_of_RNA_polII_on_the_exon_sides_of_exon_intron_and_intron_exon_junctions_/372514 (Figure A and B, the junction figures). I currently use a combination of ngsplot and deepTools for most of my metagene plots, but I have no idea how I would be able to plot something like this.

Any ideas?

Answer: I have not done something like that, but using deepTools this may work (assuming you have bigwigs for the sense and antisense):

1. get a a list of exon genomic positions, be sure that they contain the orientation.
2. For a figure similar to A, use computeMatrix reference-point, select as reference point the end of the region and 150 bp before and 150bp after the reference point.
3. For a figure similar to B do the same, except that you choose as reference point the start of the region.
Question: We are in the process of developing analysis pipelines for WGS,RNA and Methyl-seq data. The first projects considers 100 patients with a common disease. All the analysis would be run on data from these patients. Our institute's computational center is offering us 5 nodes (24 CPU core, 2.6GHz clock speed, 128GB RAM memory each), I believe that necessary storage space would also be provided to us. Based on your experiences would you please tell me if having 5 nodes is going to be sufficient, I understand that there cannot be one single answer to this question but any help from you all is much appreciated.

Answer: Since you appreciate that there can't be one right answer here you are off to a good start.

Your requirements will/may change over time and you should build some flexibility into the allocation to access nodes with better specs (mainly RAM) in case you find that you are unable to complete some specific analysis. This recommendation assumes that no *de novo* assembly work would be involved since the specs (mainly memory) would be insufficient for that type of analysis.