cannin/biostars_qa huggingface dataset:
Question: There are two VCF files that I like to merge them, using GATK or VCFtools. The problem is, they have different chromosomal notation, one has Chr, the other does not. This question could be similar to [this one](/p/13462/)

Is there any quick awk/sed commands that you suggest ?! Also I appreciate if you make comment, which of these two (GATK/VCFtools) is more reliable for this task.

Answer: Another way of doing this by using bcftools v1.9:

echo "1 chr1" >> chr_name_conv.txt
echo "2 chr2" >> chr_name_conv.txt
bcftools annotate --rename-chrs chr_name_conv.txt original.vcf.gz | bgzip > rename.vcf.gz
Question: I know for the reverse procedure, there are lots of methods available. But how can I change my fasta file with sequences like:

>gi|189094002
MDDAELNAIRQARLAELQRNAAGGGSSTNPSSSSSGGAQDSAQENMTITILNRVLTNEARERLSRQQTKITFNRKNIAADDDEDDDDFFD
>gi|68485955
MDDAELNAIRQARLAELQRNAAGGGSSTNPSSSSSGGAQDSAQENMTITILNRVLTNEARERLSRVKIVRRDSQQKQQTKITFNRKNIAGDDEDDDDFFD

to like:

>gi|189094002
MDDAELNAIRQARLAELQRNAAGGGSSTNPSSSSSGGAQDSAQENMTITILNRVLTNEARERLSR
QQTKITFNRKNIAADDDEDDDDFFD
>gi|68485955
MDDAELNAIRQARLAELQRNAAGGGSSTNPSSSSSGGAQDSAQENMTITILNRVLTNEARERLSR
VKIVRRDSQQKQQTKITFNRKNIAGDDEDDDDFFD

I know its just about inserting line breaks in each sequence, but is there any way that I can check whether the fasta file is single line fasta or not, and if is, add line breaks in each sequence..

Ps. I am using windows..

Thanks for your consideration

Answer: Since the OP has asked for a solution that will work on **windows**, `reformat.sh` from BBMap can be used like so (NN = number of characters you want on one line)

reformat.sh in=your_file.fa out=new_file.fa fastawrap=NN

On windows you would use it like so (adjust windows file path for BBMap install)

java -ea -Xmx200m cp=C:\bbmap\current\ jgi.ReformatReads in=x.fa out=y.fa fastawrap=NN
Question: How to extract allel, Genotype from vcf file using python or other language for 23GB files?
Well, I am able to right script to get allel but for large VCF file its difficult ? what should other possible way to get allel, Genotype information?

Answer: See [bcftools query][1].

----------
**EDIT:** WIth `bcftools query` you can print any information you like. So in your case e.g.:

$ bcftools query -f '%CHROM %POS %REF %ALT [ %GT]\n' input.vcf

The output looks now like this:

chr1 10177 ACC ACCC 0/1
chr1 10327 T C 0/0
chr1 10352 TAC TAAC 1/1
chr1 12783 G A 1/1

fin swimmer
[1]: http://www.htslib.org/doc/bcftools.html#query
Question: In legacy blast there is an option `-F`, which turns on or off filtering. By default filtering is ON and it effects my results a lot, so I turn it off by `-F F` (false)

In blast+ I couldn't figure out the equivalent of `-F F`. In another post it was suggested to use `-dust yes`, I tried it and it does not work, there's still filtering going on.

So, what is the option in blast+ to turn off filtering?

Answer: Seconds later I posted the question I noticed my mistake, I was supposed to test `-dust no`.

When I use `-dust no` filtering is off.
Question: I launched samtools to call SNPs on 2 BAM files, and when it comes to BCFtools, both returned an error.

$ samtools mpileup -ugf refseq.fa aln.bam | bcftools view -bvcg - > var.raw.bcf
$ bcftools view var.raw.bcf | vcfutils.pl varFilter -D 100 > var.flt.vcf

It returns the following error:

[bcf_sync] incorrect number of fields (6 != 5) at 13:51289753

What does it mean? Can I just skip this error and keep going with my list of SNPs and indels?

Answer: So I found my way around the problem. When samtools hasn't finished running the output BCF files are truncated which triggered that error when I run BCFtools.
Question: I need a way to assign a number to genes representing how broken it is depending on how many SNPs with minor alleles appear on the gene. Based on just a handful of SNPs though, not every SNP.

So lets say for example I know a few SNPs for COMT:

```
Rs165599 AA
Rs6269 AA
Rs737865 AA
Rs4633 TT
Rs769224 GG
Rs4680 AA
```

and some SNPs for STAT3:

```
Rs744166 GG
Rs8069645 AG
Rs9891119 CC
Rs6503691 CC
```

and I wanna find out which gene is more defective (as in which has the most diminshed ability to encode the corresponding protein), how can I do that? These genes obviously have way more SNPs than just those ones, and in every gene I will have many minor alleles so I'm guessing its only certain SNPs that cause problems for the gene. So could I find a number that tells me how much a bad allele of a particular SNP breaks the gene and then use that number to get an overall idea of how malfunctioning the gene is?

Answer: There is a number of suggestions on this thread.
Question: Hi All,

Is it possible to re-run STAR RNA aligner to obtain wiggle plots on already aligned and sorted BAM files? Or do I need to re-align again to obtain this type of output? I really do not want to have to re-align - so I am open to other options.

Thank you!

Answer: `bamCoverage` from deepTools would be the quickest method to do this. You can also have it output bigWig rather than wiggle files, which is generally preferable (you can directly load those in IGV, for example).
Question: I am using the lastest MutSigCV v1.4, and trying to implement the ["CoMut plot"][1] to visualize the results. However, the [MutSigCv documenttation][2] only describes the content of one the results files (`my_results.sig_genes.txt`), but there are 5 files that MutSigCV outputs:

- `my_results.categs.txt`
- `my_results.coverage.txt`
- `my_results.mutations.txt`
- `my_results.mutcateg_discovery.txt`
- `my_results.sig_genes.txt`

Yes, I did ask the developers but... *tumble-weed*

**Q1: Where can I find an explanation for the data in the other 4 files**?

As I mentioned, I am replicating some parts of [coMut][1], which has generated a couple of questions in Biostar ([here][3] and [here][4]). The plots representing the top mutated genes (left of coMut) and type of mutation per sample (bottom) were easy using data from MutSigCV - using data in the files `my_experiment.sig_genes.txt` and `my_experiment.mutations.txt`

However without having a precise explanation of the MutSigCV output, I cannot proceed further.

**Q2: Which data is used for the remaining components of the plot**?

I have a suspicion that the top bar plot (number of mutations/Mb) uses data from [MuTect][5], but no idea about the other 2 components. My guess is that some of that is buried in the MutSigCV output files but can't figure out where.

[1]: http://www.broadinstitute.org/blog/visualizing-cancer-genome
[2]: http://www.broadinstitute.org/cancer/software/genepattern/modules/docs/MutSigCV/1?print=true
[3]: https://www.biostars.org/p/13226/
[4]: https://www.biostars.org/p/87350/
[5]: http://www.broadinstitute.org/cancer/cga/mutect

Answer: There is now an R package that reproduces this type of plot, and much more:

[GenVisR](https://github.com/griffithlab/GenVisR)

May thanks to the Griffith Lab. Other options listed [here](https://www.biostars.org/p/77286/).
Question: Hi there,

I am currently working on a project related to pan-cancer analysis. I have done differential expression analysis for miRNAs and genes with edgeR. I know that edgeR only takes raw counts as input, so I downloaded the HTSeq-Count data from GDC data portal. After the differential expression analysis, I'd like to obtain the normalized expression values of miRNA/genes (RPKM/FPKM) for the downstream analysis, such as using pearson correlation between miRNAs and mRNAs to construct a miRNA-mRNA regulatory networks and so on. However, I got stuck for days on how to get the normalized expression values of both miRNAs and genes. Here are my questions:

1. There is a function called "cpm"(counts per million) in edgeR, but it says it doesn't take the gene length into accout; edgeR also provides another version of normalized counts "pseodu.counts", however, someone says this is quite difficult to interpret. So I am wondering if I could use "logCPM" as the normalized expression values for the downstream analysis?

2. If not, I realized that there is also a function called "rpkm" in edgeR which could calculate the normalized expression values. However, it needs the gene/microRNA length information to make it work. I do not know where to find the length information for genes and microRNAs, since there is no such information contained in the HTSeq-Count file. Could any one please tell me how to do it? Should I download the gene information from ENSEMBLE and the miRNA information from mirbase? And calculate the length information by myself? Is there any R package that could do the work instead?

3. Could I just download the RPKM files of miRNAs and genes from GDC data portal to construct the miRNA-mRNA regulatory network? Is that right? It seems to be the easiest way for me to do though....

Any help would be really appreciated.

Answer: There is no much use of calculating RPKM/FPKM for miRNAs.

To answer your question, you can use [featureCounts][1] to quantify your genes using a GTF file. This outputs a matrix of counts, which also includes a column of gene lengths. This column can be given to edgeR `rpkm()` function.

Option 3 would make sense. If you already have normalised data available, you can use to calculate the correlations.
[1]: http://bioinf.wehi.edu.au/featureCounts/
Question: Hi all,

I need some help with understanding the output of CNVkit, specifically Segmented log2 ratios (.cns) and the exported CNVs in VCF format.

I'm looking to get the copy number of every region found by CNVkit.

For Segmented log2 ratios (.cns) file:

Please correct me if i'm wrong: To get the actual estimated copy number I should simply anti-log the log2 column, right?

In that case,

- What's the correlation/connection, if there is any, between the log2 value in the .cns file and the inferred SVLEN value in the .vcf file?

- Is there any connection to the "CN" (copy number genotype...) value in the .vcf? Also, why does "CN" only appears in duplication events? I tried calculating the copy number from the log2 value in the .cns file but the values are different than what I expected.

**An example if it helps:**

For a certain region in the .cns file, the log2 value is 0.382191.
In the .vcf file, SVLEN is 6812878 and the CN value is 3. What is the copy number then?

Thanks,
Alon

Answer: For segmented log2 ratios -- nearly. See the documentation for CNVkit's [call][1] command. But the easiest way to get the integer copy number values is through the "export" command, either VCF or BED.

- I recommend ["export bed"][1] for custom analysis. If you are analyzing tumor samples with some known amount of normal-cell contamination, use "call" or "rescale" first to adjust the log2 ratios for that contamination.
- The ["export vcf"] output conforms to the VCF spec, which is a little unintuitive for describing CNVs. The CN field means absolute copy number, but it's only specified for copy number gains, while hemizygous or homozygous deletions are specified differently, as you're seeing.

In your example, a log2 value of .38 corresponds to `2^(.38) = 1.3` times the reference ploidy. For a diploid genome, the absolute copy number would be `2 * 2^(.38) = 2.6`, which you can probably round up to 3 assuming some reasonable level of normal-cell contamination. (See the documentation on [tumor heterogeneity][4] for more guidance on this topic.)

[1]: http://cnvkit.readthedocs.org/en/latest/pipeline.html#call
[2]: http://cnvkit.readthedocs.org/en/latest/importexport.html#bed
[3]: http://cnvkit.readthedocs.org/en/latest/importexport.html#vcf
[4]: http://cnvkit.readthedocs.org/en/latest/heterogeneity.html
Question: I'm trying to assemble small (20Mb), diploid fungal genome from MiSeq reads (~400bp after merging, 100x coverage).

The tricky thing: it's heterozygous. The divergence is ~4%, but there are hundreds (thousands?) loss of heterozygosity (LOH) regions, accounting in total to almost half of the genome...

Do you know of any assembler (methodology) capable of handling with such data?

I have tried many assemblers/tricks:

- de Bruijn graph: Velvet, ABySS, SOAPdenovo
- overlap-based: MIRA, Newbler
- clustering reads and then assembling with CAP3

but with rather bad effects. Every time, the assembly was very fragmented (N50 ~10kb), homozygous regions (LOH) were collapsed with ~200x, while heterozygous regions were separated with 100x.

I will be happy to hear about any ideas:)

**EDIT**

As no satisfactory solution exists, I'm developing [Redundans][1] (check [below][2] for more info).

[1]: https://github.com/lpryszcz/redundans
[2]: https://www.biostars.org/p/63047/#139529

Answer: As there were no good solutions, I've developed a tool that assists assembly of heterozygous/hybrid genomes ([Redundans][1]).

The program takes as input **assembled contigs**, **paired-end and/or mate pairs sequencing libraries** and returns **scaffolded homozygous genome assembly**. Resulting assembly should be **less fragmented** and with total **size smaller** than the input contigs. In addition, Redundans will automatically **close the gaps **resulting from genome assembly or scaffolding ([more details][2]).

Hopefully, someone else will benefit from that! I'll be really grateful for your feedback.

BTW: sorry for being necromancer;)

[1]: https://github.com/lpryszcz/redundans
[2]: https://github.com/lpryszcz/redundans/blob/master/test#redundans-pipeline
Question: In a lot of biological oriented paper, one often finds
figure like this:

What is the program usually people used for drawing that?
If R, how do they achieve the horizontal marker and asterisks?

Answer: Not the most elegant solution, but you can use the regular R utilities for drawing lines and asterisks. You can capture the center positions of the bars by assigning the barplot call to a variable. The lines() and text() functions take x and y coordinates of where to draw things:

# capture x coordinates of bars
x <- barplot(c(1,4,3,6,5))
# create the y coordinate of the line
y <- 4.5
# set an offset for tick lengths
offset <- 0.2
# draw first horizontal line
lines(x[1:2],c(y, y))
# draw ticks
lines(x[c(1,1)],c(y, y-offset))
lines(x[c(2,2)],c(y, y-offset))
# draw asterics
text(x[1]+((x[2]-x[1])/2),y+offset,"**")

You simply have to do some guessing at the offsets.
Question: Hi all.

I am trying to understand the mutation frequency to annotate my exome sequence data.

Since I look it up on the Internet but I tried to fail how to calculate it.

To make it matter worse, I didn't know the actual meaning of mutation frequency (Wiki seems to explain their definition but I didn't get it..it is too formal) and why this information is important when analyzing mutation.

Could anyone tell me a little bit about their basic meaning and equation to get it?

Answer: Frequency is just a count. It is usually normalized to some fixed unit of time or space to enable comparison with other counts. Mutation frequency can conceivably refer to many things depending on the context:

- Number of mutations per sample/per Mb/per gene, etc
- Number of samples in which a particular mutation is observed
- Percent of reads that support a particular mutation
- Number of mutant alleles in an individual or population (usually called "allele frequency")

Really, you need to define for yourself what it is you want to calculate. The calculation itself should then be straightforward.
Question: I have a problem building tophat 2.1.0 on Debian jessie. Running configure does not seem to produce errors (see below) but running make causes fails with

```
#some output omitted here

Makefile:1347: recipe for target 'reads.o' failed
make[2]: *** [reads.o] Error 1
make[2]: Leaving directory '/home/pommy/inst/tophat-2.1.0/src'
Makefile:313: recipe for target 'all-recursive' failed
make[1]: *** [all-recursive] Error 1
make[1]: Leaving directory '/home/pommy/inst/tophat-2.1.0'
Makefile:253: recipe for target 'all' failed
make: *** [all] Error 2
```

Here are the configuration results

- tophat 2.1.0 Configuration Results

```
C++ compiler: g++ -Wall -Wno-strict-aliasing -g -gdwarf-2 -Wuninitialized -O3 -DNDEBUG -I./samtools-0.1.18 -pthread -I/home/pommy/include -I./SeqAn-1.3
Linker flags: -L./samtools-0.1.18 -L/home/pommy/lib
BOOST libraries: -lboost_thread -lboost_system
GCC version: gcc (Debian 4.9.2-10) 4.9.2
Host System type: x86_64-unknown-linux-gnu
Install prefix: /home/pommy
Install eprefix: ${prefix}
```

Answer: Solution found here: https://github.com/infphilo/tophat/pull/13

I managed compiling a fork that uses a newer SeqAn version:

```
git clone https://github.com/wookietreiber/tophat.git
cd tophat
git checkout update-seqan
./bootstrap
./configure
make
make install
```

(See also https://groups.google.com/forum/#!topic/tuxedo-tools-users/M7Z9zLuu8o4)
Question: We are hoping to use 1000 Genomes samples as a population control for our study. The 1000 Genomes Project provides fastq, BAM and VCF files on their ftp site. We do not want to use VCF files as they have been filtered and might not contain variants occurring in our samples (especially false-positive variants in our samples). Using dbSNP is problematic for the same reason.

So it seems like a good alternative is to use 1000 Genomes BAM files. However, it would save us compute time if we could use gVCF files. Does anyone know if gVCF files from 1000 Genomes Project samples are publicly available?

Answer: A few years later, gVCF files from the 1000 Genomes project are now publicly available on [Figshare][1]

[1]: https://figshare.com/collections/1000_Genomes_gVCFs/4414307/3
Question: I just do some practice on ballgown and stringtie, and I got some GTF file and ballgown`s file.

However, I just find I can't use R or something else to deliver a excel file which contain FPKM.

I think the excel file I want maybe looks like this:

gene_id FPKM
A 124
B 541
C 122

Please help me

Thanks a lot :)

Answer: Form StringTie output you can use 'sample1_**gene_abund.tab**' file to extract these information,

**FPKM:**

cat {sample1}_gene_abund.tab | cut -d$'\t' -f1,8 | sed "s/FPKM/sample1_FPKM/" | sed 's/Gene ID/gene_id/' > {sample1}_gene_abund.fpkm.txt

**TPM:**

cat {sample1}_gene_abund.tab | cut -d$'\t' -f1,9 | sed "s/TPM/sample1_TPM/" | sed 's/Gene ID/gene_id/' > {sample1}_gene_abund.tpm.txt
Question: Hello,

I would like to get the coverage for every single base in all regions in a bed file. For this I use samtools now

`samtools depth -aa -b regions.bed input.bam`

This all works as excpected. But the more regions I have the longer this tasks takes of course. I can speed it up by dividing the bed file and start parallel processes. But it still needs about 45min to complete. My whole analyse pipeline from fastq to an annotated vcf file runs only 1h for this data set.

So I wonder if there is a faster way for this task?

fin swimmer

Answer: Try bedtools coverage with -d option ([link][1])

(And do notice the first note on *-sorted*)

[1]: http://bedtools.readthedocs.io/en/latest/content/tools/coverage.html
Question: Hello clever community!

I need your advice.
I am working on a *de novo* plant genome assembly of ~400 Mb.
I have Chromium 10x data, which was assembled with supernova. I also have Illumina paired end reads. Now I have additional data of PacBio reads, 120x roughly. The genome is diploid and I am thinking about using Falcon.

What do you think should be the best strategy:

1. Assembling PacBio reads and then using a tool to integrate the two assemblies? Is there anything like this? Which tool would you use?

2. Using a tool that can assemble the genome from both the chromium and the PacBio reads? Is there anything like it?

3. Assembling the PacBio reads and using chromium 10x and the illumina for polishing? If I assemble with Falcon, what tool should I use for polishing?

4? Anything else that I am missing to get the best out of what I can get?

Thank you very much in advance!
Alex

Answer: I don't have any experience with Falcon, so I can't help on that, but this is my advice on your other questions. In any case, your PacBio coverage is quite decent, so you might expect relatively good results from the Falcon assembly.

1. Integration of assemblies is usually not trivial - though I might just lack a good reference.
2. [SPAdes][1] can integrate all that data afaik. You can use the chromium contigs as "*trusted contigs*"
3. No experience, but my guess is you risk to polish out any heterzygosity.
4. [MaSuRCa][2] can handle both PacBio and Illumina - however, you won't be able to use the chromium data directly.

[1]: http://cab.spbu.ru/software/spades/
[2]: http://www.genome.umd.edu/masurca.html
Question: Hello

I have imputed data from ukbiobank in bgen format. I would like to convert it to a vcf file.

I can use plink2 to make pgen files and then use plink2 again to create a vcf

plink2 --bgen ukb_imp_chr17_v3.bgen --sample ukimp_chr17_v3_s.sample --make-pgen

plink2 --pgen plink2.pgen --pvar plink2.pvar --psam plink2.psam --export vcf

This creates a vcf file but it doesn't seem to process in any of our pipelines.

1. Does what I did look correct?
2. How do I check the vcf file for accuracy?

Thank you in advance ,

Richard

Answer: Assuming you want dosage information in your VCF, you need to replace "--export vcf" with something like "--export vcf vcf-dosage=DS". You may also want to add the 'bgz' modifier to request bgzipping of the VCF file.
Question: Hello

I have a .ped file in plink and would like to convert it into the 012 format (0 hom ancestral), 1 het, 2 dom derived. Is there any way to achieve this with an option in plink?

Answer: plink --file [.ped/.map fileset prefix] --recodeA --out [output prefix]
Question: Hello,

I downloaded mutation data from TCGA from here: http://firebrowse.org/ select SKCM >SNP6_copynum

and the data looks as follows:

Sample	Chromosome	Start	End	Num_Probes	Segment_Mean
TCGA-3N-A9WB-10A-01D-A38I-01	1	61735	534551	31	-0.3442
TCGA-3N-A9WB-10A-01D-A38I-01	1	564621	16345255	8476	0.0159
TCGA-3N-A9WB-10A-01D-A38I-01	1	16363940	16389566	8	0.8209

Is there a way I can find out which type of mutations these are? I tried looking at the SNP6 manifest and could not get that information either.

Thank you!!

Answer: It is processed (by GISTIC 2.0, I believe) copy number data derived from Affymetrix SNP 6.0 microarray, which has probes that target both SNPs and copy number variants. It is a commonly used platform for determining genome-wide CNV profiles (I used it during my PhD in breast cancer).

Note that you may want to instead download the file with '*minus_germline_cnv*' in its name, i.e., download it from Firebrowse.

If you want to process it, take a look here:

1. Part I
2. Part II
3. Part III

Kevin
Question: Please help.

I have two files (file1 and file2). I would like to extract the columns from file2 that have their IDs listed in file1. These are big files, with thousands of columns and lines.

file1

Id123B

Id124A

Id125A

file2

Code sex id123B id127 id125A

desired output file:

id123B id125A

Code i have tried:

#!/usr/bin/perl
use warnings;
use strict;
use feature qw{ say };

open my $COLUMNS, '<', shift or die $!;
chomp( my @columns = <$COLUMNS> );

open my $DATA, '<', shift or die $!;
my @header = split /\t/, <$DATA>;
my %column_index;
@column_index{ @header } = 0 .. $#header;

@columns = grep exists $column_index{$_}, @columns;

while (<$DATA>) {
chomp( my @cells = split /\t/ );
say join "\t", @cells[ @column_index{ @columns } ];
}

Answer: Your code should work fine, but you are not including the headers in the output. Just add `say join "\t", @columns;` before the while-loop.

BTW, `split /\t/` is equivalent to `split`. :)