cannin/biostars_qa huggingface dataset:
Question: I am using HOMER's mergePeaks tool (homer/v4.6) to find the overlapping peaks between different bed files, for downstream visualization. However, the resulting number of peaks do not match the number of peaks that I started with. I have files that look like this:

# a set of ChIP-Seq peaks
$ wc -l Sample1.bed
106536 Sample1.bed

$ head Sample1.bed
#chr	start	end	name	conc
chr1	10009	10438	chr1:10009-10438	1.93238645557258
chr1	710223	710919	chr1:710223-710919	2.6191877990312
chr1	712387	715259	chr1:712387-715259	6.36598230427597
chr1	752136	752850	chr1:752136-752850	2.77465732982152
chr1	755212	755710	chr1:755212-755710	2.15283141685661
chr1	756766	759022	chr1:756766-759022	4.05288454015793
chr1	760940	763718	chr1:760940-763718	6.44867373331963
chr1	772008	781353	chr1:772008-781353	6.86725737719796
chr1	800093	801813	chr1:800093-801813	4.41658071582566

# TSS regions from Gencode
$ wc -l gencode.bed
105785 gencode.bed

$ head gencode.bed
chr1	1868	22010
chr1	19553	39554
chr1	20266	40366
chr1	42472	62473
chr1	43048	63049
chr1	52947	72948
chr1	59090	79091
chr1	121024	141025
chr1	150445	170446
chr1	307719	327730

And the command I am using is this:

mergePeaks Sample1.bed gencode.bed -prefix mergepeaks -venn venn.txt

The stdout stream even displays the correct number of peaks:

mergePeaks Sample1.bed gencode.bed -prefix mergepeaks -venn venn.txt
Max distance to merge: direct overlap required (-d given)
Merging peaks...
Comparing Sample1.bed (106535 total) and Sample1.bed (106535 total)
Comparing Sample1.bed (106535 total) and gencode.bed (105785 total)
Comparing gencode.bed (105785 total) and Sample1.bed (106535 total)
Comparing gencode.bed (105785 total) and gencode.bed (105785 total)

However, the resulting files do not contain the full number of peaks.

$ cat venn.txt | cut -f3-
Total	Name
13229	gencode.bed
40836	Sample1.bed
19649	Sample1.bed|gencode.bed

13229 + 40836 + 19649 = 73714 peaks, which is less than the 100k that both sets started with. This is also reflected in the line counts for these files:

$ wc -l mergepeaks_*
40837 mergepeaks_Sample1.bed
19650 mergepeaks_Sample1.bed_gencode.bed
13230 mergepeaks_gencode.bed
73717 total

Any idea what is happening to these missing peaks?

Answer: After a thorough investigation, I ran mergePeaks on the original gencode.bed file by itself and found that it contained a number of peaks that overlapped each other, thus reducing the total number true peaks by the same amount as were missing. This was corroborated by running `bedtools merge -i` on the file and producing output with the same number of entries. Mystery solved. The original HOMER outputs were indeed accurate.
Question: Hello,

I am trying extract certain information from a gbk file I can extract the locus tag and the amino acid sequence however I am struggling to extract the gene location as it not in the same format in the file e.g.: `/locus_tag="NCTC86_00002"`

This is my script so far:

```py
from Bio import GenBank
from Bio import SeqIO

gbk_filename = "HS.gb"
faa_filename = "HS_converted.faa"

input_handle = open(gbk_filename, "r")
output_handle = open(faa_filename, "w")

for seq_record in SeqIO.parse(input_handle, "genbank"):
print "Dealing with GenBank record %s" % seq_record.id
for seq_feature in seq_record.features :
if seq_feature.type=="CDS" :
assert len(seq_feature.qualifiers['translation'])==1
output_handle.write(">%s from %s\n%s\n" % (
seq_feature.qualifiers['locus_tag'][0],
seq_record.name,
seq_feature.qualifiers['translation'][0]))

output_handle.close()
input_handle.close()
print "Done"
```

Answer: Hi,

the SeqFeature objects have a "location" attribute that contains the start/stop position of the feature.

```py
from Bio import SeqIO

gbk_filename = "HS.gb"
faa_filename = "HS_converted.faa"

output_handle = open(faa_filename, "w")

for seq_record in SeqIO.parse(gbk_filename, "genbank") :
print "Dealing with GenBank record %s" % seq_record.id
for seq_feature in seq_record.features :
if seq_feature.type=="CDS":
assert len(seq_feature.qualifiers['translation'])==1
output_handle.write(">%s from %s\n%s\n" % (
seq_feature.qualifiers['locus_tag'][0],
seq_record.name,
seq_feature.qualifiers['translation'][0]))
print('Start: %d, Stop: %d, Strand: %d'%(int(seq_feature.location.start),
int(seq_feature.location.end),
seq_feature.strand))

output_handle.close()
print "Done"
```

Hope this helps
Question: How to analyse Level 3 data (RPKM values) from TCGA? What Kind of analysis could be performed using level 3 data to find out differential expressed genes between normal and diseased samples?

Answer: For RNA-seq data you can use DESeq2 or EdgeR to perform a differential expression analysis. These tools are part of Bioconductor in R. However, both of these programs perform their own internal normalizations and they recommend you input the raw counts, not the RPKM or RSEM scaled estimates (for RNAseqV2 TCGA data). If you are downloading from the TCGA Portal each patient's data is in a separate file, but if you go to https://confluence.broadinstitute.org/display/GDAC/Dashboard-Stddata and click on the Open link beside the cancer of interest you will find tar files that contain merged text files with all patients in one file. Once you get all the raw counts or normalized counts in a single matrix you can analyze the data with any program that accepts a matrix of data; BUT make sure it is meant to be used on RNA-seq data because this type of data has different properties than microarray data and needs to be treated slightly differently when genes have low or zero read counts. In doing my own comparisons between DESeq2 and EdgeR, I have found I prefer DESeq2 results because it compensates for low reads counts, which can artificially inflate fold changes.
Question: Hi all,

I have multiple files that have the same (long) title format:
Lactococcus_lactis_subsp_cremoris_strain_number_file.gbk_fasta.fasta_proteins.fasta_results

I want to delete everything after ".fasta" so it would be something like:
Lactococcus_lactis_subsp_cremoris_strain_number_file.gbk_fasta.fasta

I tried to make a small bash script in order to do that quickly but when I run it, nothing is changing.
Here is my script:

#!/bin/bash

for file in /home/Documents/Folder;
do
cut -d_ -f1,2,3,4,5,6,7,8 <<< "$f";
done

There is obviously something wrong with it but I can't find where ...

If you have any ideas/suggestions for improvements, please, let me know.

Thank you for your help!

Answer: `cut` operates on file contents, not file names - not unless used much differently. Plus, your loop variable is called `$file` while your loop uses something called `$f`.

You should go with sed. It could be done with bash parameter expansion too, but sed is easier.

In your case, you'd need something like

echo "mv ${file} $(echo $file | sed 's/.fasta.*/.fasta/')" #untested - the first "." in the sed expression might need escaping

If the above `echo`s `mv` commands as expected, remove the echo and the double quotes surrounding the `mv` command to execute the command.
Question: Does anyone know how to download and install the R-GSEA package?

If I try to register via http://www.broadinstitute.org/cancer/software/gsea/software/software_index.html, I just get a 404 error.

Is there any other source for this package?

Answer: There are several packages implementing GSEA on Bioconductor, gage and phenoTest are two examples, surely someone will point more.
Question: Hi,

I am using RMarkdown in Rstudio and I want to execute commands from a program using a bash chunk

````{bash}
```

I have a program called samtools on my computer so when I execute it in the chunk, it works.

```{bash}
samtools
```

When I type :

```{bash}
which samtools
```

The output is that it tells me samtools is located in the usr/local/bin directory.

However, when I execute a program with vcftools , I get an error because Rstudio, does not know where the program is:
I have it in another directory on my computer.

How do I get Rstudio or Rmarkdown to execute vcftools from the bash chunk?
Is there a way that I can tell RMarkdown which directory to look in to find the program?

For example (something like) :

```{bash}
$vcftools = /Users/m.o.l.s/Programs_For_Bioinformatics/vcftools

```

or would I have to move all of the programs to usr/bin/local?

Outside of Rstudio, I have made aliases to the programs so they work fine on the terminal.

I made the alias by writing in my .bash_profile

`alias bcftools=/Users/paths/to/where/the/program/is/installed`

but I added the path to vcftools to my export PATH in the bash profile completely.

Answer: `.bash_profile` is executed when you start a login shell. My suspicion is that the `bash` R markdown cells are executed in a separate shell instance, that is not a login shell and so this is not executing your `.bash_profile` file when it starts.

It might be different with `.bashrc`, which might get executed when you start a bash shell in Rmarkdown (or it might not - for example, im pretty sure neither are run on an SGE submission script).

You could try explicitly adding `source ~/.bash_profile` to the start of your code chunk.
Question: Hi.
I am new for performing RNA-Seq. I have got the sequencing raw reads (PE) for my samples from Arabidopsis thaliana plants. I have analysed these reads using FASTQC. Now I have to do the alignment and mapping of these reads with the refrence genome.
In this regard my first question is, should I do any further cleaning or processing for my reads before alignment. And should I have to merge the pair-end reads?
Second question is which aligner should I use? I have tried with BWA, Bowtie2 and STAR, and I got the maximum alignment with STAR.
Other question is which source I should use for the reference genome. Like TAIR10, Araprot11 or ftp://ftp.ensemblgenomes.org/pub/release-41/plants/fasta/arabidopsis_thaliana/dna/

I will be grateful if you may guide me for my analysis.

Best
Umesh

Answer: You should :

- Not do further cleaning of the fastq files
- Not merge paired-end fastq files
- Use STAR for alignment (BWA and bowtie2 are maint for DNA not RNA)
- For the genome you can use the ENSEMBL version : https://plants.ensembl.org/Arabidopsis_thaliana/Info/Index
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: As suggested by Pierre, BLAST will pretty much answer all of your questions. You can use the LinkOut Links (usually on the right, use CTRL+F to find them) on the NCBI webpage to get additional information.

You should note one potential pitfall though: in order to properly search for proteins, your contigs should come from mRNA sequencing, because eucaryotic DNA is likely to contain introns in their sequence (which you can't make sense of that easily). If you're searching genomic sequences, finding the right BLAST parameters can be a bit tricky.

In most cases you will know a bit of information about your sequence that will help you applying the right searches. If you're going totally blind, doing a BLAST against NCBI's nr database is a good place to start (because most likely your sequence is not entirely new or unknown).

What you can try as well:

1. Trying to find an ORF; mRNAs proteins will have a possible translation from nearly the beginning to nearly the end; Polyadenylation signals (in eucaryotes) might provide additional help (there are 3 STOP codons- think about how likely it is for a given DNA sequence of length N to contain none by pure chance)
2. Use NCBI's Conserved Domain Search to classify the protein
3. There are many specialized databases to find information as soon as you know what type your sequence is; for instance, there is BRENDA for enzymes or RDP for rRNA.
4. If it is a protein and you are interested in structure, you might want to take a look at the PDB.

Sorry if this was too basic, but you mentioned you were new to biology ;-)
Question: Hi everyone, i have a blast p result with the average identity percentage 35%, is it an acceptable percentage?

If not what is the minimal acceptable identity percentage?

Between Identity percentage, e-value and bitscore, which one should we focus on in order to find the best match? Thanks.

Answer: Identity 35% means that 35% of aa in your sequence match to other sequences in database. There isn't something like "acceptable percentage". It always depends on what you are looking for:
--- if you have unknown protein sequence and you would like to know the homology sequences, information about identity (even 35%) is valuable,
--- if you have known protein and you need to confirm the sequence, the identity 35% is small and may suggest that something went wrong during your analysis.

The E-value is very important, the lower the better.

Best,

Agata
Question: Hello! Is it possible to retrieve information about variants of a single gene from the ExAC Browser in VCF format? I can only see the option to "Export table to CSV".

Answer: You can use tabix and query directly from the FTP:

tabix -h ftp://ftp.broadinstitute.org/pub/ExAC_release/current/ExAC.r0.3.1.sites.vep.vcf.gz 2:39967768-39967768 > exported.vcf
Question: Hi - I am exploring how many threads I can use for psiblast. The default value is 1 thread, and I want to know how high I can set this. Using

$ cat /proc/sys/kernel/threads-max
126335

I obtain the maximum number of processes available to the system. I somehow doubt that this is analogous to the amount of threads I can use in the psiblast `num_threads` flag. can anyone shed any light on this?

Thanks

Answer: Don't go near 126 thousand; the max threads for an operating system includes moving the hard-drives and networking, and everything else. See "ps aux" for a list of how many things are already active!

Choosing a threading level for a program is based on a dozen factors, usually if you want to do it if you have more CPU cores sitting idle, while the single-threaded version of psiblast is running 100% of one core. Try two and see if you can get two cores to 100%. Maximally you'd benefit from the number of threads = number of CPU cores, unless there is other non-cpu constraints.

It depends on the system architecture and the software architecture. In my experience, a dual processor, quad core (for 8 threads max) can do something like compute pi in eight threads efficiently; If instead of computing pi, you were reading big files from disk, and have only the one disk, youll see any more than 1-2 threads slowing each other down as they have to wait for data. Running six threads accessing the harddrive will be slower than three. This ratio depends on how much of each resource is needed by each thread.

Something like sequence alignment (BWA or Bowtie etc) needs to read a little data, then crunch a lot of CPU, so spinning up all 8 threads is fine, theyll wait their turn for data, and then get off-sync from each other and end up with 100% disk utilization and almost 8x100% CPU.

If your processes are for example requesting data from a web-server with some unknown delays, then you could run a dozen or a hundred threads and theyll wait and go when they can.

It also depends on the motherboard architecture. You probably only have 2 or 4 channels to access RAM, so running more than 4 threads that need high volume access to RAM will also slow each other down.

The key word is contention. The usual solution is trial and error; you will have to measure the speed for various threads-settings and choose an optimal for your task. I dont know about psiblast specifically, but it probably needs to access RAM quickly, and youll see it get slower per thread after 5 threads. Maybe the optimal is 6-8, but I guarantee trying 1,000 simultaneously will not be faster than 10.

Finally, of course is the level of parallelism available to the algorithm, sometimes BLAST has to work sequentially and will ignore your thread setting for some parts of the job, so these optimal settings can vary with the reference genome and query sequences.
Question: I am looking for a tool/script/pony to correct the REF column in a vcf file whenever that nucleotide doesn't match the reference genome, as supplied in fasta format.

It sounds like a common task but I could not find something. I did find `bcftools +fixref` but that only works for SNPs. My vcf files are from structural variants.

Cheers,
W

Answer: Hello https://www.biostars.org/u/24526/ ,

there is another `bcftools` plugin that should do this:

$ bcftools +fill-from-fasta input.vcf -- -c REF -f genome.fa > output.vcf

fin swimmer
Question: Hi Everyone,

I am working on Whole Genome Sequencing and analysis of Human genome from illumina HiSeq platform with about 30X coverage. Each sample (human genome) have about 250-300 fastq.gz files, whom I am dealing with 'fastqc' for quality check using following command :

/usr/local/bin/fastqc -t 8 -f fastq -o OUT/ -casava *.gz -noextract

Although it is running fine and generating equal number of "fastqc.zip" files which I unzipped using unzip '*.zip'. So, here I have 2 questions:

1. Can I merge two or more fastq files and then run fastqc on those merged files? If yes, how should I merge those fastq files?
2. I have to manually check 250-300 fastqc folder to know the quality by opening .html page. Is there any way where I can have summary of overall quality of the fastq files in a flowcell?

Please let me know your comments. I'll be highly thankful to you.

Best,
Ravi

Answer: We have a script that will run fastqc and generate a [summary report][1] with the images from all the fastq files it was run on. You may also find it useful to systematically parse the fastqc_data.txt files from each run and combine the results that way.

The script is here, but may not be the most useful and documented thing ever... Depends on imagemagick to generate thumbnails...

https://github.com/metalhelix/illuminati/blob/cluster/scripts/fastqc.pl

Also uses this script:

https://github.com/metalhelix/illuminati/blob/cluster/scripts/thumbs.sh

[1]: http://research.stowers.org/mcm/fastqc/fastqc_plots.html
Question: Hi guys! I need help in setting a `for` loop in R (I'm quite new programming in R). I would like to add the same header to all the files that match a concrete pattern inside a folder.

To get the list of files I'm using the following code:

filelist <- list.files(pattern = "DESeq2_result*")

And this is the `for` loop I am trying to implement:

for (i in seq_along(filelist)) {
names[[i]] <- a
out [i]
}

where `a` is a vector that I defined with the names of the different columns:

a <- c("gene_id", "baseMean", "log2FC",
"SD", "WaldStatistic", "pval", "padj")

If you have any tutorial/page to help me to learn and practise my ability to code functions and loops in R would be so much appreciated.

Thank you so much in advance!

Jordi

Answer: Use `col.names = ` argument when reading the files, then write out, something like this:

for(i in list.files(pattern = "DESeq2_result_*"))
write.table(read.table(i, col.names = c("gene_id", "baseMean", "log2FC",
"SD", "WaldStatistic", "pval", "padj")), i)

**Note:** this overwrites existing files, to create new files:

for(i in list.files(pattern = "DESeq2_result_*"))
write.table(read.table(i, col.names = c("gene_id", "baseMean", "log2FC",
"SD", "WaldStatistic", "pval", "padj")),
paste0(i, ".fixed.txt"))
Question: Trying to index vcf file but getting the following

```
tabix -p vcf dbsnp_138.hg19.vcf.gz
Not a BGZF file: dbsnp_138.hg19.vcf.gz
tbx_index_build failed: dbsnp_138.hg19.vcf.gz
```

Thoughts on how to proceed? Thanks!

Answer: Looks to me like the dbsnp file is not bgzipped?

```
gunzip dbsnp_138.hg19.vcf.gz
bgzip dbsnp_138.hg19.vcf
tabix -p vcf dbsnp_138.hg19.vcf.gz
```
Question: How can I get the number of mapped reads for a particular region?

`samtools view -c -F 4 my.bam` gives me count in the entire bam file but I can't just add `-r Chr1:0:1000` to get reads in that region only.

Answer: samtools view in.bam chr1:0-1000 | wc -l

but better way is to do it with bedtools if you have many regions, for e.g like exon coordinates or peak coordinates, which is more efficient way for counting purpose.

bedtools multicov -bams aln1.bam [ aln2.bam aln3.bam . . ] -bed ivls-of-interest.bed
Question: Dear all:

I have downloaded Geneious software trial version [https://www.geneious.com/][1], but now it is expired. As a fresh level grad student, I can't afford to pay full version of Geneious software. I aware that Geneious has abundant features to analyze NGS data. My expectation is to find tools that relatively similar to Geneious. My major intent is to carry out **phylogenetic analysis**.Can anyone recommend the better alternative for Geneious tools? Is there any free, open source software tools that can be alternative for Geneious? My colleague is not well comfortable with R environment, so we are seeking tools that require less programming input. Can anyone give possible aid or point out which library we can use?

[1]: https://www.geneious.com/

Answer: Have a look at [Jalview][1]

[1]: http://www.jalview.org/
Question: Hi everyone,

I have to create an artificial .vcf file with SNVs, but also with small insertions, deletions and duplications, that will be used by tools like Exomiser, Extasy, ... for variant prioritization

For SNV's, it's quite easy to fill in the chromosome, genomic position, reference and alternative allele. However for indels & duplication, I'm often not sure how to do this correctly.

If there is a deletion for example, what is then the position that you should put? Position that is deleted, or the base before or after? Same problem for an insertion and a duplication. And which reference allele should you take in those cases? Does anyone know if there are rules to do this?

Answer: Check the VCF format specifications: http://samtools.github.io/hts-specs/VCFv4.2.pdf

In particular, in section 5 contains many examples.