cannin/biostars_qa huggingface dataset:
Question: I wrote a perl script as a wrapper for a couple of calls to binaries. One of which is a call to the EMBOSS merger as such:

$mergeCMD = "merger -asequence temp1.fa -bsequence temp2.fa -outfile temp.me -outseq temp.fa;
system($mergeCMD);

It insists on printing "Merge two overlapping sequences" for every call which annoys me to no end and clutters up my terminal. There are no evident cmd line options to turn this off, so i was thinking it could be possible to redirect this output away from stdout? Better yet would of course be if someone knew an option to make merger quiet.

Answer: You can add `-auto` to your command to turn off prompts. That is an EMBOSS standard.

The above is the right solution, but for reference, there are correct ways to capture the output of commands with Perl. [Capture::Tiny](https://metacpan.org/pod/Capture::Tiny) is lightweight does a nice job and [IPC::System::Simple](https://metacpan.org/pod/IPC::System::Simple) is a drop-in replacement for system that works better than the built-in (and avoids the shell). On that note, it is better to give a list to [system](http://perldoc.perl.org/functions/system.html) to avoid the shell and it is good to test if the command succeeded.

You can test a command by looking at the return value, but that is unfortunately not as simple as it might seem (e.g., you have to localize errors and test if an error is coming from Perl or the external command). Those are the reasons I mention these modules. They do the right thing without a lot of boilerplate or cryptic code. If you want to localize and catch the exceptions you can use [Try::Tiny](https://metacpan.org/pod/Try::Tiny). All that is just for your consideration to make life a little easier. The command will work fine as written but it is good to keep these things in mind if you are building a pipeline. It is worth mentioning that there are wrappers for EMBOSS in BioPerl, so you can do something like the following (untested).

use strict;
use warnings;
use Bio::Factory::EMBOSS;

my $usage = "perl $0 a.seq b.seq";
my $aseq = shift or die $usage;
my $bseq = shift or die $usage;

my $factory = Bio::Factory::EMBOSS->new;
my $merger = $factory->program('merger');

$merger->run({-asequence => $bseq,
-bsequence => $aseq,
-outseq => 'out.seq',
-outfile => 'out.merged'});
Question: I have some rna-seq data and they all have this odd extra peak in GC content. Adapter trimming was carried out using Skewer to produce trimmed reads. The adapter contamination is now gone but the GC profile still looks strange. Is this something I need to be concerned about? Is there a way to pull out or examine the reads that constitute this extra peak? I also noticed that the kmer profile has changed after trimming especially in the 3' end.

![fastqc plots][1]

[1]: http://s32.postimg.org/ec1w140g5/fastqc_small.jpg

Answer: In your final plots we can see that there is a sequence of DNA common to a very large proportion of your reads that occurs at the start and/or end of your reads. Trying to glue those kmers together in my head I think you end up with 2 distinct adapter sequences: "GCCGGCATCACCGCG", and "GTCGGCATCACCGCG" which are obviously both highly GC rich, but one is 2% higher than the other.

I believe this would be the cause of your other peak(s), and if you can find a way to trim that stuff off your reads, you'll be fine.
You can test this if you like with [SeQC][1], but you'd need to have a .stat method for identifying reads as starting with GCCGGCATCACCGCG. I could make one for you very quickly if you don't know python.

[1]: https://github.com/JohnLonginotto/SeQC
Question: Hi everyone :)

i want to convert my XMFA file (output of progressiveMauve) into MAF.

Has anyone had any luck converting xmfa into maf? And how did you do it?

Thanks :)

L

Answer: Thank you Seja Modha for your help :)
i found this c++ script to convert xmfa to maf:

https://github.com/Wyss/mauve-py/blob/master/mauve/src/mauveAligner/src/xmfa2maf.cpp
Question: I have fastq file with undetermined reads from Illumina HiSeq. The reads in fastq file contains illumina barcode. I do now know what barcode has been used in the lab and there is no way i could find it out. Is it possible to find it out from the fastq file I received. Are there any tools that show most repeated sub-sequence in each reads of fastq file ?

[14-01-2014] Edit1:
I have ilumina base call files (raw data from sequencer). I have possibility to convert .bcl files to .fastq using bcl2fastq provided by Illumina.

Answer: Assuming the demultiplexing was carried out specifying the index sequence location, then this could be used:

grep ^@HISEQ lane1_Undetermined_L001_R1_001.fastq | cut -d":" -f10 | uniq -c | sort -nr > barcodes.txt

This outputs a text file containing every barcode in the FASTQ file, with a count of how many times it is seen. The most frequent is at the top of the list. Could help you find the barcode?
Question: Hi all,
how can I counting SNPs and INDELs per each chromosome in a vcf file?

Answer: Using BEDOPS `vcf2bed` and standard Unix tools like `wc` to count lines:

$ vcf2bed --snvs < variants.vcf | wc -l
$ vcf2bed --insertions < variants.vcf | wc -l
$ vcf2bed --deletions < variants.vcf | wc -l

Cf. http://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/vcf2bed.html

The above does counts over all chromosomes. You need to do an extra step to count per chromosome.

To perform this count exercise quickly per-chromosome, per-variant class:

$ vcf2bed --snvs < variants.vcf > variants.snvs.bed
$ for chr in `bedextract --list-chr variants.snvs.bed`; do echo $chr; bedextract $chr variants.snvs.bed | wc -l; done

The `bedextract` tool is part of BEDOPS. If you install BEDOPS to get `vcf2bed`, you get `bedextract`, as well.

Repeat this procedure for `--insertions` and `--deletions`, replacing `snvs` accordingly, to count those classes.
Question: I am trying to perform a GSEA for ten biochemical properties to check which properties are enriched for my set of genes. I am using the Windows application (v4.2.2) with GSEAPreranked package. However, when I perform the analysis some of them returns a p-value of "0.0". When I checked the user guide it states:

"In the GSEA report, a p value of zero (0.0) indicates an actual p value of less than 1/number-of-permutations. For example, if the analysis performed 100 permutations, a reported p value of 0.0 indicates an actual p value of less than 0.01. For a more accurate p value, increase the number of permutations performed by the analysis. Typically, you will want to perform 1000 permutations (phenotype or gene_set). (If you attempt to perform significantly more than 1000 permutations, GSEA may run out of memory.)"

I tried increasing the permutation to 10,000 and 100,000 which still yielded a "0.0" value. Is there a way to calculate the absolute p-value which I can use to later calculate FDR q-value?

Answer: Hello Mahasish,

The answer to your question is implicit in the question itself.

GSEA is a form of permutation based testing. The p-value is calculated empirically by permuting the gene labels at random, then seeing how many times out of the total number of permutations a random permutation is more extreme than the actual result you generated.

example: if you run 10,000 perms, and the actual data has enrichment for a pathway more extreme than all but 13 of the permutations, then the p-value would be 13/10000 = 0.0013.

If you are generating a p-value of 0.0, this means that in 0 out of 10,000 cases was the random permutation more enriched than your actual data for that pathway.

There are two solutions I'd consider.

1) Record the p-value as p < 0.0001
2) Re-run the pathway analysis with 10,000,000 permutations. This will likely give you at least a few perms that are more enriched than your data.

If you choose 2), then I do agree with rpolicastro's comment, above: fgsea is the way to go.

Finally, are you only testing one pathway? If not, you will need to control for multiple testing.

VAL
Question: Hi everyone,

Could you please tell me where I can find the annotation and fasta file related to ERCC 92 spike ins. The Thermo fisher website only gives the fasta file. Is there any other link that I can download compatible fasta and gff files?

Thanks a lot

Answer: There is no GFF file for the spike-ins, if you need one then just add one line per spike-in to whatever GFF/GTF file you're using.
Question: Hello,

I've combined 3 seurat files in order to do an integrated analysis and everything is working fine, but I would like to extract some cells for subsequent analysis. As the dimplot is now, the clusters consist of cells from all of my conditions so I can not seem to do the usual approach which would be something like :

Cluster.1.8.12.only <- subset(x = seuratA, idents = c('1','8', '12'))
if for example I wanted clusters 1,8 and 12.

Does anyone know how I would modify the above to for example, get clusters 1, 8 and 12 but of only the condition I'm interested in? for simplicity sake lets assume the conditions are labelled A, B and C and are in an integrated seurat called ABC.combined.

Thank you for your help :)

Answer: I'd probably just subset again. Same approach, just using the metadata column using your condition.

Cluster.1.8.12.cond <- subset(x = Cluster.1.8.12.only, subset = your_condition_column == "condition_of_interest")
Question: Hi dears all,
I am going to do SNP calling using bcftools mpileup for WGS data on fish. Before SNP calling I did this preprocessing as follow:

samtools sort -n F0-F06-10.bam F0-F06-10ns.bam

# Add ms and MC tags for markdup to use later

samtools fixmate -m F0-F06-10ns.bam F0-F06-10-fixmate.bam

# Markdup needs position order

samtools sort -o F0-F06-10PS.bam F0-F06-10-fixmate.bam

# Finally mark duplicates

samtools markdup F0-F06-10PS.bam F0-F06-10-markdup.bam

Then using samtools faidx the fasta reference file got indexed and bcftools mpileup was used for SNP calling. For the last step (bcftools mpileup) I got involved with an error.

samtools faidx /home2/omid/alignment/ref/GCF_001858045.2_O_niloticus_UMD_NMBU_genomic.fa

bcftools mpileup -Ou -f GCF_001858045.2_O_niloticus_UMD_NMBU_genomic.fa.fai F0-F06-10-markdup.bam

> error: [E::fai_build_core] Format error, unexpected "N" at line 1

I will be so grateful if any one can help me to solve this problem.

Cheers

Answer: I'm quite sure that the reference file you are using in your `bcftools` command is the index file of your reference and not the reference file itself.

The index file must be located in the same location like the reference file. The file you specify with the `-f` parameter must be the reference genome file.

fin swimmer
Question: let's say I have a fasta of a protein sequence

> albumin

MKWVTFISLL FLFSSAYSRG ... ... ...

I want to split the sequence into all possible consecutive 8 amino acids (only in 1 direction, amino -> carboxyl) (And no looping(I don't know if it is the right expression), say, GMKWVTFIS)

I need

> fasta.albumin1

MKWVTFIS

>fasta.albumin2

KWVTFISL

> fasta.albumin3

WVTFISLLF

...

> fasta.albumin13

FSSAYSRG

And, I want to do this for all known human protein sequences.

How would I do it???
I need the result as a fasta or fasta files. And IDs for resulting 8-mer seuqeunces need to be unique.

Answer: There are many methods. For long terms, some programming skill is worth learning.

Here's a simple way by using [seqkit](http://bioinf.shenwei.me/seqkit/usage/#sliding).

$ seqkit sliding -s 1 -W 8 seqs.fa
>albumin_sliding:1-8
MKWVTFIS
>albumin_sliding:2-9
KWVTFISL
...
>albumin_sliding:13-20
FSSAYSRG
Question: A common error when running genemark on a fresh Ubuntu install is for gm_es.pl to fail with an error like:

Can't exec "/opt/genemark/gmes_bp_linux_64/gmes/make_mat.v3.0.ini.csh": No such file or directory at /opt/genemark/current/gmes/gm_es.pl line 270.

The paths may be slightly different, but the error will be the same. Why is this error occurring and what can be done to fix it?

Answer: The error is most likely because the c shell (csh) is not installed on your system. If this is the case and you are on Ubuntu, try:

sudo apt-get install csh

If you have csh installed in a location other than `/bin/csh`, you will need to change the first line of the `make_mat.v3.0.ini.csh` script to point to your installation, or make a symbolic link to your installation like so:

sudo ln -s /path/to/csh /bin/csh

Another cause might be that the `make_mat.v3.0.ini.csh` script is no longer executable. If that is the case, try

chmod +x /path/to/make_mat.v3.0.ini.csh
Question: I am trying to align RNASeq data to both human and mouse genome but seems getting into some trouble. I merged human genome.fa with mouse genome.fa into a combined genome.fa and then only used human genes.gtf for transcriptome assembly. Previously I have successfully aligned these RNASeq data to human genome only but my client said these RNA may have been contaminated with mouse RNA. Given there were more than 50% reads were unable mapped to human, I would try to align the data to combined human and mouse genome. Any suggestions how to do it using tophat2? With my procedure, I got the below errors:

[2016-04-13 14:06:17] Beginning TopHat run (v2.0.9)
-----------------------------------------------
[2016-04-13 14:06:17] Checking for Bowtie
Bowtie version:	2.1.0.0
[2016-04-13 14:06:17] Checking for Samtools
Samtools version:	0.1.19.0
[2016-04-13 14:06:17] Checking for Bowtie index files (genome)..
[2016-04-13 14:06:17] Checking for reference FASTA file
[2016-04-13 14:06:17] Generating SAM header for /data1/workspace/DCI/Sarantopoulos/RNASeq/Data/human_and_mouse_genome_combined/human_mouse_combined/Bowtie2Index/genome
format:	fastq
quality scale:	phred33 (default)
[2016-04-13 14:06:22] Reading known junctions from GTF file
[2016-04-13 14:06:27] Pre-filtering multi-mapped left reads
[2016-04-13 14:06:27] Mapping JP01_S9_L003_R1_001-trimmed_1 to genome genome with Bowtie2
[2016-04-13 14:41:46] Pre-filtering multi-mapped right reads
[2016-04-13 14:41:46] Mapping JP01_S9_L003_R2_001-trimmed_2 to genome genome with Bowtie2
[2016-04-13 15:16:00] Preparing reads
left reads: min. length=114, max. length=126, 32220590 kept reads (1960224 discarded)
right reads: min. length=107, max. length=126, 32456688 kept reads (1724126 discarded)
[2016-04-13 15:40:07] Building transcriptome data files..
[FAILED]
Error: gtf_to_fasta returned an error.

But there was no error when I align to human genome only previously.

Answer: I just figured out how to do it if it is necessary.

1. rename the chromosome names in mouse .fa (such as from chr1 to mchr1) so that the chromosome names will be unique across human .fa and mouse .fa
2. rename chromosome names in mouse .gtf file accordingly.
3. bowtie2 build index of mouse using renamed mouse .fa
4. merge mouse renamed .gtf with human .gtf
5. put rebuilt mouse bowtie index and human index into one directory
6. run bowtie2 for alignment and assembly. using the combined index and merged .gtf.
Question: Hello everyone,
I need plot, which visualize percentage of nucletides per base in R. I have a percentage of bases but now I don't know how to visualize them.

fastq <- readDNAStringSet("https://d28rh4a8wq0iu5.cloudfront.net/ads1/data/ERR037900_1.first1000.fastq","fastq")
freq <- alphabetFrequency(fastq)
perc <- freq/width(fastq)*100
perc <- as.data.frame(perc)
perc

Any help? Thanks.

Answer: library(Biostrings)
fastq <- readDNAStringSet("test.fa","fasta")
fastq
af=alphabetFrequency(fastq, as.prob = T,baseOnly=T)
barplot(af)

No need to convert per and as.prob=T will convert it percentages and you can bar plot directly without converting it data frame.
Question: Hello everyone,

I am trying to find a good estimate of the average and preferably also the upper and lower bounds of the mutation rate in the brown Norway rat, Rattus norvegicus. This is for simulations of rat demographic history. I have found recombination rates from a paper a few years back that compares human rat and mouse recombination rates (if anyone knows of a more uptodate estimate of rat recombination rates that is welcome too!).

Essentially while I can find a lot of work on mutation rate estimates in humans I don't see anything for the rat. If anyone can point me in the right direction that would be great. Thanks!

Best,

Rubal

Answer: This paper (http://www.sciencemag.org/content/314/5797/268), especially the supplement (http://www.sciencemag.org/content/314/5797/268/suppl/DC1) talks about calculating it. Though its human data, guess it should be the same way to apply it for rat!
Question: I have FPKM values for a number of genes. I want to make a heatmap and find differentialy expressed genes. Here I have three questions:

1. Should I use FPKM values directly to make a heatmap or is there anything which has to be done to FPKM values before using them to make a heatmap? I am using MeV to make a heatmap.
2. How do I find differentially expressed genes here? I read that taking log2 values of FPKM and then setting the cutoff as 2 will give differentially expressed gene. Taking the cutoff as 2 sounds very arbitrary, is there any reason to take 2 as the cutoff?
3. There are biological replicates also, each replicate has a specific FPKM value under each treatment. How do I treat these replicates?

Thanks

Answer: Dear Ashish,

For your situation, I would not use MeV.

I recommend reading your data into R and then transforming to Z scale via *zFPKM*. You may then use standard statistical tests on these Z-scores for the purpose of identifying statistically significant genes.

You may take a look at this thread where I give some information on producing heatmaps in R: https://www.biostars.org/p/274509/

Also, here is a more advanced version of a heatmap in R: https://www.biostars.org/p/273155/

Kevin
Question: Hi!

I have ~163000 contigs (a single FASTA file) and I want to map them to a refernce genome (also a FASTA file). Is there any way to do it? I have tried bowtie2, but could not work out its nuances.

Any help appreciated!

Thanks

Answer: Bowtie2 is a short read mapper and, although it can be used to map long sequences, it probably won't be good at it, specially if the query and reference are even moderately divergent.

In addition to the suggestions by https://www.biostars.org/u/33094/ and https://www.biostars.org/u/18713/ , there is also [LAST][1].

Another option to consider is [QUAST][2], an assembly evaluation tool. Quast uses minimap2 to align one (or more) query genomes to a reference genome, and in addition to the alignment, it will provide a number of metrics comparing the query to the reference.

[1]: http://last.cbrc.jp/
[2]: http://quast.sourceforge.net/quast
Question: Hi,

I have a couple of contigs, and I would like to answer questions of the following sort

a.) Are they unique?

b.) Which organism do they come from?

c.) Do they come from a protein encoding gene?

d.) What is that protein's function etc...

I am very new to Bioinformatics (even Bio for that matter) so any help would be appreciated. Even links to further readings.
Thanks.

Answer: The question is a bit vague, so I'll have to make a load of assumptions.

If your contigs are from a genome sequencing project, perhaps you have available raw reads as well? If so, mapping these reads might help you to determine whether your contigs are from a unique region (coverage near expected level), from a repeated region (very high coverage), or from some random contamination (low coverage). This is not absolute, of course, and you have to be careful with stuff like duplicate clones, which can artificially inflate coverage, especially if you are sloppy in the lab. If you have multiple independent libraries, something that only exists in a few of them is likely a contamination.

Another thing I've looked for in our project, is contigs that appear to be circular (having reads that map partially to both ends). The jury is still out on the usefulness of this (my one glaring example appears to be a misassembled contig).

AT/GC content might also be useful for identifying contamination, depending on the length of your contigs and the ratio in your organism.

16s and 18s sequences can classify contigs as pro- or eukaryote.
Question: There is a well-known paper in the autoimmune disease literature that describes an algorithm called probabilistically identified causative SNPs, or PICS.

The paper identified ~9,000 PICS based on a posterior probability that the SNP is a true causal variant rather than a neutrally associated variant. I trimmed this to 1,865 LD independent PICS.

I wanted to see if the PICS were enriched for having p-values below a certain threshold in our lab's GWAS data. I could have set this up as an exact binomial test using the threshold as the probability of success, but this is problematic for a lot of reasons. For instance, although the genomic inflation in our study is well-controlled (lambda GC = 1.03), that still indicates a slight enrichment of lower p-values.

For this reason and many others, I think such an exact binomial test is a weak comparison and is likely to have a very high type I error rate. So, I searched for a more rigorous comparison.

Without going into substantial detail, after a lot of work I think a more appropriate test is to compare the association p-values of the PICS (from our study) to p-values of "SNPs like them" (e.g., SNPs that have similar numbers of nearby genes, or similar numbers of exonic variants, etc. etc.).

So, using SNPsnap I generated a list of 87,000 SNPs having characteristics matching those of the PICS. I call them "matched SNPs."

I then drew random samples of 1865 from from the ~87000 in the matched SNP list and compared the number of times the matched SNPs had more variants with an association p-value below the threshold than the PICS. I did this 100,000 times, and each time the matched SNP sample had more SNPs below the threshold than the PICS SNPs, I recorded this.

So, if 4650 of the 100,000 samples from the matched set had more SNPs below the p-value threshold than the PICS set, then I am saying p = 0.04650 and calling that an "empirical p-value."

I do not call this a permutation based test anywhere, because I have not scrambled any case-control statuses or anything like that. However I am thinking of calling it a "simulation test." Also, I would like to ask if it is valid to call this an "empirical p-value," and if not, I am unsure what to call it.

My question is, does either of these terms carry with it a denotation that would make it inaccurate to use that term in the above context?

Answer: Your general strategy sounds reasonable and yes this is an empirical p-value (do mention the simulation in the methods!). The appropriateness of the result is dependent on how "similar" the simulated data was, so I hope you put a good deal of thought into that. Regarding the p-value itself, that's the correct value assuming none of the simulations returned the exact same number of SNPs below threshold, since those should be included too.

**Edit:** Apparently there's an [argument to be made](http://www.ncbi.nlm.nih.gov/pmc/articles/PMC379178/) that 1 should be added to the numerator and denominator when calculating the p-value. I've never thought about it, but it seems reasonable at first read.
Question: Trying to get site that have a FILTER value of "PASS" or "."

I know of:

`bcftools view -f PASS filename.vcf.gz`

and

`bcftools view -f . filename.vcf.gz`

I tried:

`bcftools view -f 'PASS | .' filename.vcf.gz`

but no dice.

Is it possible to get what I am looking for or do I have to filter separately then concat?

Answer: bcftools view -f 'PASS,.' filename.vcf.gz