cannin/biostars_qa huggingface dataset:
Question: Hello biostars!

I was wondering, how to compare a few samples in my multisample vcf file?

Basically, I would like to get a venn diagram, but vcf-compare is working only on multiple vcf files, bcftools stats and plot-vcfstats do not help either (they do not compare positions across the samples).

I can just split my vcf on three separate files and compare them using the approach described above, but I believe there should be more simple solution.

Answer: Using [Bioalcidae][1] and the following script?

var counts={};
while(iter.hasNext()) {
var ctx = iter.next();
for(var i=0;i< ctx.getNSamples();++i)
{
var gi = ctx.getGenotype(i);
for(var j=i+1;j< ctx.getNSamples();++j)
{
var gj = ctx.getGenotype(j);
if( gi.sameGenotype(gj) ) {
var key = gi.getSampleName()+"-"+gj.getSampleName();
var c= counts[key];
if(c==null) c=0;
++c;
counts[key]=c;
}
}
}
}

for(var i in counts) {
out.println("["+i+"]\t"+counts[i]);
}

run:

gunzip -c input.vcf.gz| java -jar dist-2.0.1/bioalcidae.jar -F vcf -f script.js | column -t

[S1-S10] 35
[S1-S2] 16
[S1-S3] 20
[S10-S2] 16
[S10-S3] 20
[S2-S3] 19
[S1-S4] 16
[S10-S4] 16
[S3-S4] 19
[S2-S4] 15

[1]: https://github.com/lindenb/jvarkit/wiki/BioAlcidae
Question: Hello,

I am not sure if this has been answered before as I looked and couldn't find a simple answer.

I have a bam file, and all I want is to annotated all regions with 0 coverage in bed format. Is that possible?

Thank you,

Adrian

Answer: You want the ["genomecov" tool](http://bedtools.readthedocs.org/en/latest/content/tools/genomecov.html#bga-reporting-genome-coverage-for-all-positions-in-bedgraph-format)

Example:

$ bedtools genomecov -ibam aln.bam -bga | head
chr1 0 554304 0
chr1 554304 554309 5
chr1 554309 554313 6
chr1 554313 554314 1
chr1 554314 554315 0
chr1 554315 554316 6
chr1 554316 554317 5
chr1 554317 554318 1
chr1 554318 554319 2
chr1 554319 554321 6

To extract solely the intervals with no coverage just use awk to filter the results:

$ bedtools genomecov -ibam aln.bam -bga | awk '$4==0' | head -n 2
chr1 0 554304 0
chr1 554314 554315 0
Question: Hi guys,

As shown below, I have two files a fasta file (containing a sequences and sequence IDs) and a text file containing a list of sequence ID.

I would like to extract the respective sequence for each of the sequence IDs in the text file. I have tried a number of times, but I simply cannot even open the fasta file - I'm a total beginner with regards to bioruby - if you could explain any answers in detail, I would be highly grateful.

EDIT: Could you give your answers using Ruby or BioRuby. Thanks

e.g. EXAMPLE FASTA FILE

>CATH_RAT
MWTALPLLCAGAWLLSAGATAELTVNAIEKFHFTSWMKQHQKTYSSREYSHRLQVFANNWRKIQAHNQRN
HTFKMGLNQFSDMSFAEIKHKYLWSEPQNCSATKSNYLRGTGPYPSSMDWRKKGNVVSPVKNQGACGSCW
>CATH_MOUSE
TFSTTGALESAVAIASGKMMTLAEQQLVDCAQNFNNHGCQGGLPSQAFEYILYNKGIMGEDSYPYIGKNG
QCKFNPEKAVAFVKNVVNITLNDEAAMVEAVALYNPVSFAFEVTEDFMMYKSGVYSSNSCHKTPDKVNHA
>CATH_HUMAN
MNPTLILAAFCLGIASATLTFDHSLEAQWTKWKAMHNRLYGMNEEGWRRAVWEKNMKMIELHNQEYREGK
HSFTMAMNAFGDMTSEEFRQVMNGFQNRKPRKGKVFQEPLFYEAPRSVDWREKGYVTPVKNQGQCGSCWA
>CATH_MONKEY
FSATGALEGQMFRKTGRLISLSEQNLVDCSGPQGNEGCNGGLMDYAFQYVQDNGGLDSEESYPYEATEES
CKYNPKYSVANDTGFVDIPKQEKALMKAVATVGPISVAIDAGHESFLFYKEGIYFEPDCSSEDMDHGVLV

EXAMPLE .TXT FILE

>CATH_MOUSE
>CATH_HUMAN

Many Thanks

Answer: This is what I ended up doing - by using 'pure' ruby.

```
require 'bio'
# create a hash of the Fasta file
def analyse_input_file(input_file)
input_read = {}
biofastafile = Bio::FlatFile.open(Bio::FastaFormat, input_file)
biofastafile.each_entry do |entry|
input_read[entry.entry_id] = entry.seq
end
input_read
end

# Iterate ID file and look up in hash
def iterate_id_file(id_file, input_read, output_file)
ids = File.read(id_file)
ids.each_line do |id|
# remove the inital > and the end of line character.
id = id.gsub!(/^>/, '').chomp
seq = input_read[id]
File.open(output_file, 'a+') do |f|
f.puts ">#{id}"
f.puts seq
end
end
end

# Set up
input_file = ARGV[0] # '/Volumes/Data/input_file.fa'
id_file = ARGV[1] # '/Volumes/Data/ids.txt'
output_file = "#{input_file}.out"

# Set up
input_read = analyse_input_file(input_file)
iterate_id_file(id_file, input_read, output_file)
```

Just a quick explanation on how the above works... (I know I would have found this really helpful, when I asked the question)

The main concept used here is the HASH (also known as a dictionary).

The basic arrangement of a Hash is as follows

Hash_name = { 'key_1' => 'value_1', 'key_2' => 'value_2' }

This is also written as:

Hash_name = { key_1: 'value_1', key_2: 'value_2'}

The main benefit of Hashes is that you can quickly call a value by passing the respective key as follows

```
puts Hash_name[key_1]
# => value_1

puts Hash_name[key_2]
# => value_2
```

In the `analyse_input_file` method, we created a hash where the key was the sequence id (i.e. `entry.entry_id`) and the sequence itself was the value (i.e. `entry.seq`)

Then in the `iterate_id_file`, we simple passed the id to the HASH to get the corresponding sequence....

Simples...
Question: Hi all,

I have calculated the `MEDIAN_CV_COVERAGE`, however, what the meaning of this measurement? I have checked the [related discussion at seqanswers][1] and asked the question there. In brief, if the `MEDIAN_CV_COVERAGE` as the defined on wiki, the calculated value should be larger than the output of picard. In my data sets, this value just range between 0.2 and 1.7. To my understand, this value will definitely larger than 0 for normal samples. Because the coverage for high expressed genes should be very high and the variation for the top 1000 highly expressed genes will be large. Am I right?

[1]: http://seqanswers.com/forums/showthread.php?t=21420&highlight=MEDIAN_CV_COVERAGE

Thank you!

Best,
p

Answer: The `MEDIAN_CV_COVERAGE` represents a measurement of the variability of coverage over single genes. Ideally, there would be no variation over the length of the genes, but in reality, there is always some. Extreme outlier samples may represent technical differences in those samples that lead to higher variability. In practice, take a look in a browser (IGV, etc.) at the extreme samples on the low and high end to get a sense of what the value represents.
Question: Here's another beginner BioPython question from me...

I'm running some genome assemblies for someone who has some new Illumina sequence data and also had done some sequencing a few years ago. They have some Sanger and 454 sequences (a couple thousand sequences with a couple thousand base pairs for each) and everything is in fasta format. They are trusted contigs and I'd like to use the data to help bridge assembly gaps of Illumina sequencing data (from the same extracted DNA that was kept frozen for a few years -- same bacterial strain). There are no quality files for the older data and if they ever existed they are long gone from the data directories.

What I want to do is create FASTQ files for the old Sanger and 454 FASTA files with a phred score of 40 -- I think that would be a `I` under PHRED-33.

(I know there are assemblers out there that take fasta files inputs along with fastq -- but my in-house makefile and pipeline does not and I really would like to incorporate these sequences in our current workflow without a lot of changes).

There are lots of resources for converting from fasta to fastq *WITH* a qual file -- see [here](/p/16458/), [here](/p/85929/), [here](http://qiime.org/scripts/convert_fastaqual_fastq.html), [here](http://seqanswers.com/forums/showthread.php?t=16925), and [here](http://sequenceconversion.bugaco.com/converter/biology/sequences/fasta_to_fastq-solexa.php); but I've looked all over but can't find anything for converting without a qual file.

[Here is a public perl script](https://code.google.com/p/fasta-to-fastq/) that does exactly what I want it to do -- but I can't get it to work for me and I keep getting an error.

As an exercise, I quickly wrote the script below, but I can't get it to provide the quality score as the length of the sequence. The quality string needs to be the same length as the input sequence string. Anyone have any hints for me? Any help is appreciated.

Here's my code so far:

#!/usr/bin/env python

"""
Convert FASTA to FASTQ file with a static quality read

Usage:
$ ./fasta_to_fastq NAME.fasta NAME.fastq
"""

# import libraries
import sys, os
from Bio import SeqIO
from Bio.SeqIO.QualityIO import PairedFastaQualIterator

# Get inputs
fa_path = sys.argv[1]
fq_path = sys.argv[2]
phred_quality = "I"

# Check inputs
if not os.path.exists(fa_path): raise Exception("No file at %s." % fa_path)

# make fastq
with open(fq_path, "w") as handle:
records = PairedFastaQualIterator(open(fa_path), print("phred_quality"))
count = SeqIO.write(records, handle, "fastq")

# print fastq
print "Converted %i records" % count

Answer: The solution is to create a record within each loop and assign qualities to that. This has some explanations: http://biopython.org/DIST/docs/api/Bio.SeqIO.QualityIO-module.html

```
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio import SeqIO

seq = Seq("NACGTACGTA")
rec = SeqRecord(seq, id="foo")
rec.letter_annotations["phred_quality"] = [40] * len(seq)
recs = [ rec ]
handle = file("out.fq", "wt")
SeqIO.write(recs, handle, "fastq")
```

Better yet do it with bioawk:

```
$ more fasta.awk
{
qual = ""
for (i=0; I
Question: Is anyone aware of a tool which allows one to apply the concept of "mutation hotspots" to 3D protein structures? It is common to define mutation hot spots at the DNA level where we define the hotspot as a region of DNA sequence with recurrent mutations. But what if we want to know if some mutations are "geographically" co-located in the 3D protein structure (e.g., in the same pocket) but not necessarily "sequentially" co-located. Even without pre-defined concepts of "pockets" if a protein structure is known then it should be possible to define some concept of physical distance between residues. Any suggestions?

Answer: It looks like maybe the tool 'MuPIT' was created to help researchers address this question. From the Pubmed Abtract:

> MuPIT interactive: webserver for mapping variant positions to annotated, interactive 3D structures.

> Mutation position imaging toolbox (MuPIT) interactive is a browser-based application for single-nucleotide variants (SNVs), which automatically maps the genomic coordinates of SNVs onto the coordinates of available three-dimensional (3D) protein structures. The application is designed for interactive browser-based visualization of the putative functional relevance of SNVs by biologists who are not necessarily experts either in bioinformatics or protein structure. Users may submit batches of several thousand SNVs and review all protein structures that cover the SNVs, including available functional annotations such as binding sites, mutagenesis experiments, and common polymorphisms. Multiple SNVs may be mapped onto each structure, enabling 3D visualization of SNV clusters and their relationship to functionally annotated positions. We illustrate the utility of MuPIT interactive in rationalizing the impact of selected polymorphisms in the PharmGKB database, somatic mutations identified in the Cancer Genome Atlas study of invasive breast carcinomas, and rare variants identified in the exome sequencing project.
Question: Dear all,

considering a RNA-seq experiment and analysis that provides the expression values as TPM, please would you let me know what is a minimum TPM value in order to consider a gene to be expressed ?

talking about RPKM.FPKM units, I remember that a gene was considered expressed if RPKM (or FPKM) > 1 ... thanks a lot,

-- bogdan

Answer: There is no such thing as a cut off, because there is no such thing as not expressed - the whole genome is transcribed at some level in any given cell type. However, this doesn't stop us from sometimes needing to make a decision: which genes to include in an metagene analysis for example.

One a purely technical note, one could define a cut off as the point at which we can't distinguish between low expression levels and technical noise. It seems like zFPKM is doing something similar to this, but hand-rolled versions I've seen shift exons into nearby, but unexpressed, GC matched genome regions and then quantify them to get an average signal for strictly unexpressed sequence. I guess the difference between this and zFPKM is that zFKPM gives you a level for "unexpressed genes", where as this gives out a baseline level for "not genes".

A different approach is to think about what the level actually means. This is next to impossible with FPKM, but FPKM is readily translatable into TPM, and then the meaning is quite concrete. For example, if the average cell has 200,000 mRNA molecules in it at any one time, then a TPM of 5 would translate to 1 molecule per cell on average at any one time.

Finally, you could think in a distributional sense. I think in the last sample I looked at TPM 5 put you in the top 10% most highly expressed genes.

In the end it depends on what the purpose of the threshold is.
Question: Hi,

I'm looking for a nice way to calculate the percentage of reads on specifique chromosomes (MT und (1-23,X,Y) und unplaced scaffold).

I know I can get the read counts per chromosome with `samtools idxstats`, however I have a lot of bam files and I would like to automate the calculation. The problem is, I'm struggling with basic batch text maniplulation and would therefore appreciate any help you can give me (specific or general).

Edit: I forgot to mention, that I have multiple bam files, which are indexed and only consist of uniquely mapped reads.

Answer: General workflow you could try given that you use R.

1. Look at using `dir()` to get a list of index files.
2. Look at using `sprintf()` to create a command string to run samtools idxstats given a filename (which you have from #1), piping the results to a unique file for each sample.
3. Look at using `system()` to run the command string from 2.
4. Look at using a for loop loop to run #3 over your list of index files from #1.
5. Define a function to read a single idxstats output file.
6. Loop using sapply or lapply over all the idxstat output files
7. Combine output of 6 into a single matrix, with samples in columns and chromosomes in rows.
Question: Dear BioStars Leaders,

I have been running BWA MEM 0.7.12 on a handful of paired-end sample fastq files and it worked out pretty well. I have also included BWA as part of a bash and python based pipeline.

When I run BWA on a specific lane of a specific sample fastq files (paired-end sequencing files), I received this error and wondering what the issue might be.

[gzread] : invalid distance code

Here is the command that I ran on a Linux server :

```
$bwa mem \
-t 18 \
-M \
-R "@RG\tID:development_run_070_WES-VAL3_L002\tSM:Sample_13016\tPL:IlluminaNextSeq500\tLB:Lib1\tPU:Unit1" \
/home/hg19/ucsc.hg19.fasta S4_L002_R1_001.fastq.gz S4_L002_R2_001.fastq.gz > bwaAlignReads.sam 2> bwa.stderr.log
```

I am curious to hear from others if you have got similar error, and it would be great if anyone could suggest any possible solutions

Here are the last 30 rows from the bwa-mem error log file that I saved :

```
[M::process] read 1314136 sequences (180000260 bp)...
[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (2, 531348, 12, 3)
[M::mem_pestat] skip orientation FF as there are not enough pairs
[M::mem_pestat] analyzing insert size distribution for orientation FR...
[M::mem_pestat] (25, 50, 75) percentile: (134, 177, 227)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 413)
[M::mem_pestat] mean and std.dev: (182.75, 71.19)
[M::mem_pestat] low and high boundaries for proper pairs: (1, 506)
[M::mem_pestat] analyzing insert size distribution for orientation RF...
[M::mem_pestat] (25, 50, 75) percentile: (125, 227, 409)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 977)
[M::mem_pestat] mean and std.dev: (210.27, 167.43)
[M::mem_pestat] low and high boundaries for proper pairs: (1, 1261)
[M::mem_pestat] skip orientation RR as there are not enough pairs
[M::mem_pestat] skip orientation RF
[M::mem_process_seqs] Processed 1309276 reads in 506.988 CPU sec, 28.235 real sec
[W::bseq_read] the 2nd file has fewer sequences.
[M::process] read 774742 sequences (106186594 bp)...
[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (2, 533411, 8, 4)
[M::mem_pestat] skip orientation FF as there are not enough pairs
[M::mem_pestat] analyzing insert size distribution for orientation FR...
[M::mem_pestat] (25, 50, 75) percentile: (134, 177, 226)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 410)
[M::mem_pestat] mean and std.dev: (182.10, 71.36)
[M::mem_pestat] low and high boundaries for proper pairs: (1, 502)
[M::mem_pestat] skip orientation RF as there are not enough pairs
[M::mem_pestat] skip orientation RR as there are not enough pairs
[M::mem_process_seqs] Processed 1314136 reads in 497.150 CPU sec, 27.629 real sec
[gzread] : invalid distance code
```

Thanks

Answer: can you please run :

gzip --test S4_L002_R1_001.fastq.gz S4_L002_R2_001.fastq.gz
Question: As part of an accepted paper I would like to share both the code and data files so readers can reproduce my results. The code is easy to share via a public repository (e.g. GitHub). However, I am unsure about how to share the data. This is a computational project, so we use published genomic sequence data, process it (mapping/filtering etc), and analyze. I would like to provide the users with the processed data files prior to the downstream analysis (the downstream analysis can be run using the code I will provide). So my question is:

**Where can/should I deposit the processed data files for sharing?**

Many public scientific data repositories declare that they only accept new data which has not been published already. But the datasets I use have already been published (at least as raw data, or processed with a different pipeline).

Just to be clear, the motivation here is to save users the need to reprocess the entire raw data, which could require significant effort and computational resources.

Answer: We post the datasets we use for Galaxy trainings (fastq, bam, bigwig, vcf, etc.) to Zenodo. This is generally convenient enough, since it's free.
Question: I sorted my RNA-Seq bam files after mapping by star, and here is my key command line in a for cycle bash file:

samtools sort $file -o ${filename}.sorted.bam

At last, the one error file shows like:

[bam_sort_core] merging from 6 files and 1 in-memory blocks...
[bam_sort_core] merging from 4 files and 1 in-memory blocks...
[bam_sort_core] merging from 8 files and 1 in-memory blocks...
[bam_sort_core] merging from 9 files and 1 in-memory blocks...
[bam_sort_core] merging from 7 files and 1 in-memory blocks...
[bam_sort_core] merging from 4 files and 1 in-memory blocks...
[bam_sort_core] merging from 5 files and 1 in-memory blocks...
[bam_sort_core] merging from 3 files and 1 in-memory blocks...
[bam_sort_core] merging from 6 files and 1 in-memory blocks...
[bam_sort_core] merging from 13 files and 1 in-memory blocks...
[bam_sort_core] merging from 8 files and 1 in-memory blocks...
[bam_sort_core] merging from 2 files and 1 in-memory blocks...
...

Something wrong here? I searched but still don't know why, so I'm not sure these sorted files be used for next step?
What's more, I have M(like 63) files, and there M-N(like 6, yes, few) lines, it means not all file will meet this problem? It makes me more confused.

Any ideas will be appreciated!!!

Answer: These are not error messages, just debugging notes.

Large files cannot be sorted in memory thus get saved into temporary files. Once the sort completes the temporary files are removed.

There is nothing to be concerned about
Question: I understand that Refseq is a NCBI project. But I wonder if there is anyway to obtain the data from Ensembl?

I find Ensembl seems a little more organized at times but I'm used to NCBI and not at all familiar with Ensembl interface. I am considering fiddling around with the Ensembl interface.

Answer: See the table at the bottom of this view for RefSeqs mapped to Ensembl transcripts:

http://www.ensembl.org/Homo_sapiens/Gene/Matches?g=ENSG00000139618;r=13:32889611-32973805

Using the Perl API, you can find these as 'Xrefs' or External References.

If you want to use the browser, have a look at the videos on our youtube channel, or go through our 'website walkthrough':

http://www.youtube.com/user/EnsemblHelpdesk

http://www.ensembl.org/info/website/tutorials/coursebook.pdf

More tutorials can be found here:

http://www.ensembl.org/info/website/tutorials/index.html

Feel free to email our helpdesk, particularly for browser-related questions:

helpdesk@ensembl.org
Question: Hello everyone!

I am trying to learn more about the SAM/BAM format and how reads are mapped to the genome.

I understand that reads can be mapped to the forward (sense) and reverse (antisense) strand.

I understand that in every read pair there is a first-in-pair and a second-in-pair read - and that this probably has something to do with the order in which they came off the sequencer and not the position in the genome? (but maybe im wrong)

I understand that weird things can happen. A pair of reads can map millions of base pairs apart. Both reads can map to the same strand. But can reads map in both upstream and downstream direction - i.e. the read sequence once reversed (but not reversed complimented) was mapped.

With three variables thats 32 different combinations for paired reads (and another 4 for singletons).

With two variables, there's only 8 different combinations (and 2 for singletons)

Regardless, does anyone know of any tool available which will give you the statistics for the 32+4 or 8+2 combinations in your BAM file? If not, i will add it to the new version of metaflagstat.

Here's a picture of what i'm trying to explain :)

Answer: I was planning on posting this along with the two programs that uses these read types (a QC program which counts these read types, and a pileup tool (like bedtools, etc) where you can filter on read type - but since they're a long way from being published I'll just post the list for now :)

EDIT: I should add that the formula on the right are also not 100% true. I fell into the trap of assuming the read end is the read start plus the length of the sequence. Its not - you have to parse the CIGAR string. So wherever you see `len(SEQ)`, that really means `len(of_all_the_mapped_and_skipped_bases)`. But that didn't fit. Also the variable PSEQ doesn't exist, which means some values cannot be calculated without knowing the mate's read data too.