cannin/biostars_qa huggingface dataset:
Question: Do freeBayes and platypus support gVCF?

https://sites.google.com/site/gvcftools/home/about-gvcf

With GATK you used to have the 1000 + 1 genome variant calling problem, where you had to redo the whole expensive multi-sample variant calling (gatk haplotype caller) if you wanted to add a single genome.

Now GATK Haplotype caller can create intermediate gVCF files for each sample and do a much less expensive merging of the gVCF files to a multi-sample vcf.

https://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_gatk_tools_walkers_variantutils_GenotypeGVCFs.php

Do freeBayes and and platypus also support this or are they planning to support this?

Answer: My understanding is that freeBayes, platypus and samtools do not yet support gVCF files. However, you can work around this by calling variants from individuals files, merging the variant calls and then recalling the individual bam files at all sites called variant in at least one individual. See [this excellent blog post][1]

[1]: http://bcbio.wordpress.com/2014/10/07/joint-calling/
Question: I had 421 vcf files, I want to convert the all vcf files to the ped file (format of plink),

and I want to merge these 421 ped files to a ped file, what is the best should I do?

Thanks a lot!

Answer: using [plink][1]
in the directory that contains your vcf files run this command "assuming that plink in the [path][2]"

for i in *.vcf; do f="${i%.*}"; plink --vcf $i --maf 0.05 --recode --out $f".ped";done

example:

plink --file fA --merge-list allfiles.txt --make-bed --out mynewdata

[1]: http://pngu.mgh.harvard.edu/~purcell/plink/
[2]: http://www.linfo.org/path_env_var.html
Question: We want use different genomic analysis tools like samtools,bcftool,plink etc

Which operating system should we use?

Answer: It really doesn't matter. Linux distributions will often have samtools in their repo., so you don't need to compile it yourself (of course, that's simple). Of course if the Unix that you install is one of the BSDs then that's probably also the case for it (likewise with Mac OS X).
Question: For my RNA-Seq project, I need a very small vcf file. I just need to use it for testing so it just has to be correctly formatted with minimal information. I was looking at the 1000 genomes project website and came up with this link: http://www.1000genomes.org/ftpsearch and I specified my search to vcf files but there are just too many to look through and attempt to download. Does anyone know of a good download from the link provided or an entirely new download I could use?

Answer: If you search for trio*vcf and scroll down a little, you get to this file:

ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/pilot_data/release/2010_07/trio/snps/trio.2010_06.ychr.sites.vcf.gz

It's 50KB, so nothing that will stress out your hardware, and you can always open it in an editor and lop off lines at the bottom. It's also got the benefit that it lists information for only a few genomes, not the proverbial 1000, so the lines are narrow enough to fit on a screen. That was my motivation to search "trio".
Question: I have multiple VCF files, and I need to copy the elements of column 12 to column 13 while keeping the column names.

From this

$ cut -f12,13 sample_file.vcf | grep -v "##" | head
VD_SCLC_1	VCFM_SCLC_1
0/1:17,17:.:.:34,13:.:.:.:7,6:.:0.2766:.:13:.:47	0/1:.:.:.:.:.:.:.:.:.:.:.:.:.:.
0/1:111:.:39:.:0.3514:.:21,18:.:.:41,31:.:.:72,39:.	0/1:.:.:.:.:.:.:.:.:.:.:.:.:.:.
0/1:.:0.4118:.:.:.:28,22:50,35:.:.:19,16:85:.:.:35	0/1:.:.:.:.:.:.:.:.:.:.:.:.:.:.
0/1:.:.:.:9,12:.:21,20:.:7,13:.:0.4878:.:20:41:.	0/1:.:.:.:.:.:.:.:.:.:.:.:.:.:.
0/1:0.3654:.:10,9:.:.:18,15:.:33,19:.:.:.:52:19:.	0/1:.:.:.:.:.:.:.:.:.:.:.:.:.:.

To this

$ cut -f12,13 sample_file.vcf | grep -v "##" | head
VD_SCLC_1	VCFM_SCLC_1
0/1:17,17:.:.:34,13:.:.:.:7,6:.:0.2766:.:13:.:47	0/1:17,17:.:.:34,13:.:.:.:7,6:.:0.2766:.:13:.:47
0/1:111:.:39:.:0.3514:.:21,18:.:.:41,31:.:.:72,39:.	0/1:111:.:39:.:0.3514:.:21,18:.:.:41,31:.:.:72,39:.
0/1:.:0.4118:.:.:.:28,22:50,35:.:.:19,16:85:.:.:35	0/1:.:0.4118:.:.:.:28,22:50,35:.:.:19,16:85:.:.:35
0/1:.:.:.:9,12:.:21,20:.:7,13:.:0.4878:.:20:41:.	0/1:.:.:.:9,12:.:21,20:.:7,13:.:0.4878:.:20:41:.
0/1:0.3654:.:10,9:.:.:18,15:.:33,19:.:.:.:52:19:.	0/1:0.3654:.:10,9:.:.:18,15:.:33,19:.:.:.:52:19:.

Answer: Try this:

grep -v '^##' sample_file.vcf | awk '$13=$12' >
Question: Of course, as a bioinformatician, I am aware of many large-scale open-source bioinformatics datasets, such as

- The ENCODE consortium (www.encodeproject.org, RNA-Seq, ChIP-Seq and so on),
- The Roadmap Epigenomics consortium (www.roadmapepigenomics.org, RNA-Seq, Chip-Seq, Bilsulfite-Seq),
- The IHEC consortium (www.ihec-epigenomes.org, RNA-Seq, Chip-Seq, Bilsulfite-Seq),
- The TCGA/ICGC consortia (www.cancergenome.nih.gov, www.icgc.org, large-scale cancer data, DNA-Seq, RNA-Seq, etc.) and
- The LINCS consortium (www.lincscloud.org/l1000 , gene expression for more than a million of different perurbation experiments).

I am wondering, however, what other wonderful datasets, the are both large and open-source, are currently available. That might include things like RNA-Seq, Chip-Seq, Bisulfite-Seq, whole genome sequencing, WGAS, and many other assays (not necessarily NGS-related, though that is what I am mostly looking for).

Also things like the (neural) connectome of certain species (in any event large data) could be of interest.

There are quite some GEO datasets that at least partially fulfill these requirements, but most are simply having to few data samples in order to be interesting to me.

Your comments are greatly appreciated!

Answer: - [FANTOM][1]
- [RegulomeDB][2] (Not a large scale, but very useful functional database)
- [1000Genome Project][3]
- [GoNL][4]

[1]: http://fantom.gsc.riken.jp
[2]: http://regulome.stanford.edu/downloads
[3]: http://www.1000genomes.org
[4]: http://www.nlgenome.nl
Question: Hi all,

I've written one perl script, which compares two vcf files and write output in third file. it is running good with small files, but when i am inputting big files in Gbs, my system becomes slow and it would hang.

I'm using dbSNP vcf file as input which is around 9GB.

Please Suggest me something that will make my perl script run faster. This is my first perl script and i'm new to perl.

Any help would be appreciated !!

Thank you !!!!

```
#!/usr/bin/env perl
use strict;
use warnings;
use Data::Dumper;
use List::MoreUtils;

open (FILE, "<", "dbSNP_in.vcf") or die "failed to open: $!\n";
my @array=();
my @CHR;
my @location;
my @rs;
my @ref_n;
my @alt_n;
foreach (@array)
{
chomp;
my ($chrom, $pos, $id, $ref, $alt, $qual, $filter, $info) = split(/\t/, $_);

push @CHR, $chrom;
push @location, $pos;
push @rs, $id;
push @ref_n, $ref;
push @alt_n, $alt;
}

open (FILE1, "<", "trial_rs.vcf") or die "failed to open: $!\n";
my @array1=();

open (OUT,">trial_output.vcf");
my @columns;

foreach (@array1)
{
chomp;
@columns=split(/\t/, $_);
my $i;
for ($i=0; $i<@array; $i++)
{
if (($columns[0] eq $CHR[$i]) and ($columns[1] eq $location[$i]) and ($columns[3] eq $ref_n[$i]) and ($columns[4] eq $alt_n[$i]))
{
$columns[2]=$rs[$i];
}

}

print OUT join("\t", @columns), "\n";
}
```

Answer: Consider the following refactoring:

```
use strict;
use warnings;
use autodie;

my ( @array, %hash );

open my $FILE, "<", "dbSNP_in.vcf";

while (<$FILE>) {
$hash{"@array[0, 1, 3, 4]"} = $array[2] if @array = split /\t/, $_, 6;
}

close $FILE;

open my $OUT, ">", "trial_output.vcf";
open my $FILE1, "<", "trial_rs.vcf";

while (<$FILE1>) {
my @columns = split /\t/;

$columns[2] = $hash{"@columns[0, 1, 3, 4]"}
if exists $hash{"@columns[0, 1, 3, 4]"};

print $OUT join "\t", @columns;
}

close $OUT;
close $FILE1;
```

- `autodie` is used to trap all I/O errors.
- `chomp` is not needed within either loop.
- `split`ting only the amount needed is faster than `split`ting all.
- An interpolated array slice (`@array[0, 1, 3, 4]`, which produces a string comprised of elements 0, 1, 3, 4) is used as a hash key with a corresponding value of array element 2--for use later.
- An interpolated array slice, as a hash key, is tested for existence, and `$columns[2]` is set to that key's corresponding value if that key exists. Doing this, instead of using four `eq` statements, should speed up the file processing.

Hope this helps!
Question: Hi!

I have to perform eQTL analysis and now I have to manipuate the vcf file. I would like to know if there are some way to convert these information: `1/1=2` ; `./.=0` ; `0/1=1`)

Example:

#CHROM POS REF ALT DGRP-038 DGRP-040 DGRP-045
2L 2262 T TTC 1/1:7:0 ./.:0:0 ./.:0:0

Then in this case the following:

#CHROM POS REF ALT sample-1 sample-2 sample-3
2L 2262 T TTC 2 0 0

Are there any software to do that?

The second question is: does the meaning of these values (0 1 2) are the following?

- reference/reference (no mutation) = 1 no present in my example
- reference/alternative (1 mutation copy) = 2 and in my case (1/1:7:0 )
- alternative/alternative (both copy mutated) = 0 and in my case (./.)

Thanks in advance!!
Tommi

Answer: With PLINK 1.9,

plink --vcf [vcf filename] --allow-extra-chr --recode A

will write 0/1/2-coded genotypes to plink.raw .
Question: Hi,
I have some mouse samples. I would like to determine sex to test for sample switch or corruption from sequencing data. I have vanila freebayes called vcf files. I figured out that there is plugin to determine sex from vcf files

bcftools +vcf2sex

However it seems to require ploidy information for given genome(see for example this problem https://github.com/samtools/bcftools/issues/175). This information should be in the format

> space/tab-delimited list of CHROM,FROM,TO,SEX,PLOIDY

I am using GRCm38.82 (http://ftp.ensembl.org/pub/release-82/gtf/mus_musculus/README) for read mapping and I cannot find that ploidy information for given genome release.

I would like to ask:

1. Which is most reliable way to determine sex of the sample either from bam or vcf file?
2. What does exactly ploidy information mean? Are those coordinates of pseudoautosomal regions? And why this information is required for vcf2sex to work properly?
3. Where can I find that ploidy information for given genome release?

Answer: I figured out that for the assemblies there are informations about pseudoautosomal regions available through Genome reference consortium. In case of my genome I had to read this http://ftp.ensembl.org/pub/release-82/fasta/mus_musculus/dna/README to understand that DNA sequence is GCA_000001635.6 which is GRCm38.p4 and then from http://www.ncbi.nlm.nih.gov/projects/genome/assembly/grc/mouse/ I can figure out which are the coordinates of pseudoautosomal regions. From that I can assemble ploidy file that might be

X 1 169969759 M 1
X 170931300 171031299 M 1
Y 1 91744698 M 1
Y 1 91744698 F 0

However I figgured out the tool works better when using ploidy file of the form

X 1 169969759 M 1
Y 1 91744698 M 1
Y 1 91744698 F 0

When facing this you can also try using this repository https://github.com/pd3/bcftools and plugin called +guess-ploidy
Question: Dear community,

I have a question how to tackle the following problem. I have a clinical dataset which has already common prognosis scores, like characteristic tumor markers, tumor classification, etc.

However I developed a new score which is able to predict patients' prognosis. But how can I show that my score is superior to other already developed prognosis markers?

I have read lots of publications with ROC-analysis and evaluation of c-indeces, ... . However, there seems to be no common consensus? Has anyone an idea how to tackle my problem?

LG,
Hashirama

Answer: Is your prognosis continuous or is it a yes/no outcome? If it is yes/no, you can make a confusion matrix and then talk about specificity and sensitivity. If it is continuous, you can measure the error (for example, root mean square error (RMSE)).
Question: Dear all,

please could you advise, what is a typical range of AF (allele fractions) of somatic mutations in the cancer samples in order to consider them clonal or sub-clonal ? Thank you ...

-- bogdan

Answer: To answer your question in a simple sentence - In an ideal scenario, variants with VAF around 50% is fully clonal. Since most somatic mutations tend to be heterozygous, a clonal mutation would have an allele frequency of 50%. Anything below ~45/40% would be a sub-clonal.

But it also depends on your tumor type. Liquid tumors (such as leukaemia) are much clean with little to no contamination. Solid tumors on other hand are dirty, often contamintated with adjacent normal tissues. In that case clonal mutations tend to be clustered around VAF lower than 50%. For example, a tumor with a purity of 80% would have a clonal mutations around 40%.

Things get messier again since one need to consider copy-number status. If a mutation is located on an amplicon, gets its VAF amplified as well. A sub clonal mutation with copy number doubled would have a high allele frequency.

There are algorithms which account for CNV and tumor purity. You can take a look at ABSOLUTE from Broad which estimates tumor purity and clonality for variants. I have seen articles, where ABSOLUTE is used to estimate cell fractions. Here is an excerpt from an article,

> Tumor cellularity was determined based on the VAF and segmented copy
> number data using ABSOLUTE, in order to determine the cancer cell
> fraction (CCF) of each mutation, as was previously described by
> McGranahan et al. The clonal status was defined according to the
> confidence interval (CI) of the CCF. Mutations were classified as
> subclonal if the upper bound of their 95% CI was less than 1.

But in general, major clone would be around 50%, assuming 100% pure sample.

Again [Chris Miller][1] is the right person when it comes to tumor heterogeneity and clonality. He could correct me or add more to it.

[1]: https://www.biostars.org/u/117/
Question: I m using deseq2 for my samples such as Wild type[untreated] , and two treated condition [with vitamin D] and Retinoic acid

and my samples are in replicates ,I want to compare the differential expression between WT vs treated[Vitamin D] and WT vs Retinoic acid .so here is my code im using but Im not sure if its doing what i want to do..

condition <- factor(c(rep("WT", 2), rep("AT", 2),rep("VD",2))))`

(coldata <- data.frame(row.names=colnames(countdata), condition))`
dds <- DESeqDataSetFromMatrix(countData=countdata, colData=coldata, design=~condition)
dds

Any suggestion or help would be highly appreciated

Answer: Given your design, you'll want to ensure that `WT` is the base level for comparisons, which it's currently not.

condition <- factor(c(rep("WT", 2), rep("AT", 2),rep("VD",2))), levels=c("WT", "AT", "VD"))

Now when you use `results()` the coefficients (2 will be AT vs. WT and 3 will be VD vs WT) will be the comparisons you want.
Question: Dear All,

I am doing few histone protein's SNP analysis.

I have used SIFT and Polyphen for my initial analysis. Now I would like predict its functional impact.

Can anyone please suggest me few tools or software for the prediction of the functional

effect of nsSNPs. Thank you

Answer: This is a good question. there are many other tools to investigate the functional effect of the SNPs.

In first instance you can model the SNPs into your protein structure of interest simple using the homology modeling (using mutate_model that you can get from the website of modeller) once you obtain your structure with the point mutation you can perform three type of analysis using three bioinformatics tools:

DSSP you can assign secondary structure both wild-type and to the variant protein and than you compare if there are a differences in the secondary structure(note that DSSP do not predict the secondary structure but only assign a possible structure)

HBPLUS you can use this tool to perform the analysis of the hydrogen bond to evaluate the pattern of hydrogen bond

NACCESS you can use this tools to evaluate the accessible areas of the molecule

If you have any questions can contact me
Question: Hi,

I have followed Trinity guidelines and assembled a denovo assembly as well as used that assembly as a reference to quantify the reads using RSEM and finally followed it with Differential expression analysis using DESeq2. For Differential expression analysis, I have used RSEM.genes.results rather than the RSEM.isoform.results, since I was not sure using isoform level expression, is accurate or not compared to gene-level expression.

But the problem now is how do I select the best transcript/isoform for the Differential expressed genes because without that I cannot extract the sequence from Trinity assembly as the Trinity assembly has sequences for isoforms and not genes.

I have thought of doing this in several ways - selecting the longest isoform, clustering all the isoforms, and then select the longest isoform but I was wondering if I can get a consensus of all the isoforms for the gene of interest.

Answer: Using isoforms instead of genes is not incorrect *per se*, but it is noisier (particularly so for a *de novo* assembled transcriptome) and needs larger sample sizes and deeper sequencing per sample, so as https://www.biostars.org/u/25721/ already said, I would indeed recommend gene-level expression analysis.

There is no need to use `tximport`, as its method has been [implemented][1] in Trinity. The **RSEM.genes.results** should be identical to importing the counts with `tximport`.

What is the "best" isoform is open to debate, but I would argue the longest is not the best. I remember seeing Trinity authors recommending selecting the most expressed isoform, but I can't find the link.

>but I was wondering if I can get a consensus of all the isoforms for the gene of interest.

It seems you want a "super-transcripts" representation of the transcripts, Trinity has implemented this as well, see the [SuperTranscripts][2] wiki page.

[1]: https://github.com/trinityrnaseq/trinityrnaseq/wiki/Trinity-Transcript-Quantification
[2]: https://github.com/trinityrnaseq/trinityrnaseq/wiki/SuperTranscripts
Question: Hi all,

Any script available for doing sequence length disrtibution of fastq file?

Thanks,
Deepthi

Answer: Save it as fastq_lenght.py and run it :)

#!/usr/bin/env python
#-*- coding: UTF-8 -*-
import sys
######################################################################################
syntax = '''
------------------------------------------------------------------------------------
Usage: python fastq_lenght.py file.fastq
result: .txt file same name as input name plus "_lenghts.txt"
------------------------------------------------------------------------------------
'''
######################################################################################

if len(sys.argv) != 2:
print syntax
sys.exit()

######################################################################################

fastq_file = open(sys.argv[1], 'r')
prefix = sys.argv[1].split('.')[0]
outfile = open(prefix + '_' + 'lenghts.txt', 'w')
seq = ''
for line in fastq_file:
line = line.rstrip('\n')
if line.startswith('@'):
if seq:
outfile.write(name + '\t' + str(len(seq)) + '\n')
seq = ""
name = line
else:
seq = line
outfile.write(name + '\t' + str(len(seq)) + '\n')
fastq_file.close()
outfile.close()

print '\n' + '\t' + 'File: ' + prefix + '_' + 'lenghts.txt has been created...'
Question: I have a multigenome vcf file. Suppose the file has samples A to Z, but I want to extract the subset of samples B to G and extract a small vcf file. How can I make such subset vcf file?

Answer: I have created this bash loop to loop over files (by chromosome or any vcf file). Then using vcf-subset tool, I was able to extract the subset file. Here, sample.txt is the list of samples per line. No need to `tabix` or `bgzip` parent vcf files with this method, but is a bit slower.

for i in /path/dir/*.vcf; do
vcf-subset -c sample.txt "$i" | bgzip -c > /get/inthis/dir/output_"${i##*/}"_.vcf.gz
done
Question: Hi everyone,

In order to polish my Falcon assembly of my D. suzukii genome, from PacBio long reads, I need to align the raw reads to my genome and create a .bam or a .sam file.

Such a tool already exist at PacBio, called blasr, but is currently not working (the 5.3 version I've installed do not provide a .bam or a .sam output... Which is VERY annoying).

Does an alternative of such aligner exist ? It has to be a aligner aware of the higher error rate of PacBio reads (that's why I can't use a classic aligner).

Any ideas ? I'm truely desesperate, I absolutely need to polish my assembly before going further in my analysis.

Thanks for your help !

Roxane

Answer: bwa and BBmap (as far as I know limited to 6kb)
___
>PacBio subreads or Oxford Nanopore reads to a reference genome:
``` bwa mem -x pacbio ref.fa reads.fq > aln.sam ```

more [here](http://www.homolog.us/blogs/blog/2013/10/23/bwa-mem-blasr-look-pacbios-newly-released-human-data/)

**Update**.

[Minimap2][1]

> is a versatile sequence alignment program that aligns DNA or mRNA
> sequences against a large reference database. minimap2 is tens of times faster than mainstream long-read mappers such as BLASR, BWA-MEM, [NGMLR][2] and GMAP.

[1]: https://github.com/lh3/minimap2
[2]: https://github.com/philres/ngmlr
Question: Hi all,

I need to do some basecalling using ont-guppy/6.1.1-1, but the configuration (flowcell and kit) is not present despite this being the last version of the software. I wonder if which other configuration would you use to replace the missing:

FLO-MIN114:SQK-LSK114

Thanks in advance!

Answer: In case someone is interested. I asked Nanopore and they said the only way is by installing the new version of ont-guppy (6.4)
Question: Hi,
I have a text file of 7500 genomic positions with their chr, start, end and want to get their nucleotide sequences. Can someone point me some tool or any thoughts on how to do it?

Thank you,

Answer: I have blog posts on this http://crazyhottommy.blogspot.com/2013/04/batch-converting-coordinates-to.html
http://crazyhottommy.blogspot.com/2015/02/fetch-genomic-sequences-from-coordinates.html
Question: I'm trying to get a regex to work with rename; I've tried the approach of similar answered questions here but couldn't get the results I wanted.

The files are named as such:

SR1_S90_L001_R1_001.fastq.gz
SR1_S90_L001_R2_001.fastq.gz
Rinc_S96_L001_R1_001.fastq.gz
Rinc_S96_L001_R2_001.fastq.gz

And I would like to retain only the information prior to the first underscore and the _R1_ or _R2_ tags, like this:

SR1_R1_.fastq.gz
SR1_R2_.fastq.gz
Rinc_R1_.fastq.gz
Rinc_R2_.fastq.gz

Thanks in advance!

Answer: Quick python solution.

#!/usr/bin/env python
import os, glob

for file in glob.glob("*.fastq.gz"):
# test with print statement
print file, '\t', file.split('_')[0] + '_' + file.split('_')[3] + '_.fastq.gz'
# uncomment to rename
# os.rename(file, file.split('_')[0] + '_' + file.split('_')[3] + '_.fastq.gz')

Save as `rename_fastq.py`; run as `python rename_fastq.py` in the directory containing fastq.gz files.

Not sure why you want to keep `'_'` after the R*
Question: Hi there,

I'm creating motif weblogo's using biopython's `motif.weblogo('somefile.png')`. It works very well, but unfortunately the resolution is super tiny (see the example) and therefore not very useful. Is there any way to crank up the resolution programatically? I'd not like to use the weblogo website as I need to produce these things from an automated script.

Thanks a lot! :-)

Answer: I'm not a python dev, but I think you should use a vectorial format like SVG rather than the default bitmap format ('PNG') when calling weblogo http://biopython.org/DIST/docs/api/Bio.motifs-pysrc.html#Motif.weblogo

motif.weblogo('somefile.svg',format='SVG')

or modify another parameter using the **kwds** argument (try stack_width="large", etc... )