cannin/biostars_qa huggingface dataset:
Question: In GWAS, what it means by 'single SNP association studies only explain a small part of disease heritability'? How this explained heritability is quantified?

Answer: In statistics, when people say "17% of the variability is explained by X Y and Z", they are referring to the proportion of the variance that can be accounted for by the predictors in the statistical model.

For example, if you did a big association study on the genetics of lung cancer. You would need to include smoking as a covariate in the model. Why? Well, 1) this will prevent you from mis-attributing cancer that is actually due to smoking to the person's genetics instead. But, your question touches on a second reason. 2) When you include smoking as a covariate you can "partial out" or more colloquially "explain away" some of the variance in the data. This means there is less total variation left over in the dataset as a whole, or, "less variance left to explain". Still confused? I'll keep yammering.

Complicated phenomena like cancer are hard to predict because they have many distinct motivators. To predict cancer incidence perfectly, you would have to know about *all* of the predictors. But if you don't, your predictive model might succeed in explaining only 50% of the phenomenon (for example). A statistician might say about this, that, "the model accounts for 50% of the variance in the dataset." In the specific context in which you are speaking, you might phrase this as "only 15% of the heritability of the disease can be explained by known genetic risk factors"...

Numerically, what is going on has to do with ratios of sums of squares. Generally, you end up with the total variance of the data. That goes on the bottom as the denominator. The amount of the variance that the model explains is the numerator, which is the amount of the total sums of squares (SS) that your model explains.

So if you had

```
(explained SS) = 10000
(unexplained SS) = 40000
```

the ratio of explained to total would be 10000/40000 = 0.25, and your model would "explain" 25% of the variance.

It is not exactly correct to say this *has to* apply (only) to complex polygenic conditions, though it is likely.

As a counterexample, consider CFTR. If mutated in certain ways, the person may develop cystic fibrosis. Now, so far, we have found >1200 mutations that can lead to that phenotype or something like it... but imagine that we had only found 300 of these mutations to date.

Despite that these variants are all in the same gene, you could still enter these 300 variants into a statistical model (e.g. general linear model) as predictors, you might be able to account for 43% of the variance in the data. This would relate directly to the ratio of sums of squares mentioned earlier.

Hope this helps.
Question: I would really like to use FastQC for my project but am getting the following error message when I try to run it on my Ubuntu server 15.04

```
bio@ubuntu:~$ fastqc &

[1] 716
rafay@ubuntu:~$ Exception in thread "main" java.awt.HeadlessException:
No X11 DISPLAY variable was set, but this program performed an operation which requires it.
at java.awt.GraphicsEnvironment.checkHeadless(GraphicsEnvironment.java:207)
at java.awt.Window.(Window.java:535)
at java.awt.Frame.(Frame.java:420)
at java.awt.Frame.(Frame.java:385)
at javax.swing.JFrame.(JFrame.java:174)
at uk.ac.babraham.FastQC.FastQCApplication.(FastQCApplication.java:71)
at uk.ac.babraham.FastQC.FastQCApplication.main(FastQCApplication.java:324)
```

Whereas java is already installed

```
bio@ubuntu:~$ java -version

java version "1.7.0_79"
OpenJDK Runtime Environment (IcedTea 2.5.6) (7u79-2.5.6-0ubuntu1.15.04.1)
OpenJDK Client VM (build 24.79-b02, mixed mode, sharing)
```

Answer: It is complaining X11 is not set (maybe not even installed?). When you call fastqc without arguments, it opens the graphical interface. Try:

fastqc *.fastq &
Question: hi,

how I can create my own annotation gtf for miRNA DE analysis by mirdeep2?

thank you

Answer: Congratulation Angel for getting the rarest badge on all of Biostars! :)

![http://i.imgur.com/HuL4eGU.png][1]

[1]: http://i.imgur.com/HuL4eGU.png
Question: I have two set data, 4 samples with 2 replicates for each one from batch 1 and another 4 samples with 2 replicates from batch 2.

I would like to remove batch effects from these samples and compare different methods together.
I have done below commands but face with error:

design

samples method batch
L4_rep1 L4 L b1
L4_rep2 L4 L b1
L6_L8_rep1 L6_L8 L b1
L6_L8_rep2 L6_L8 L b1
Q5_Q7_rep1 Q5_Q7 Q b1
Q5_Q7_rep2 Q5_Q7 Q b1
Q3_rep1 Q3 Q b1
Q3_rep2 Q3 Q b1
co_40d_A co_40d co_40d b2
co_40d_B co_40d co_40d b2
co_60d_A co_60d co_60d b2
co_60d_B co_60d co_60d b2
EB_A EB EB b2
EB_B EB EB b2
H9_A H9 H9 b2
H9_B H9 H9 b2

design$=samples <- factor(design$=samples, levels = c("L4","L6_L8", "Q3", "Q5_Q7","co_40d","co_60d", "EB", "H9"))
design$method <- factor(design$method, levels = c("L", "Q", "co_40d","co_60d", "EB", "H9"))
design$batch <- factor(design$batch, levels = c("b1", "b2"))

design.matrix <- model.matrix(~0+batch+method,design)

design.matrix

batchb1 batchb2 methodQ methodco_40d methodco_60d methodEB methodH9
L4_rep1 1 0 0 0 0 0 0
L4_rep2 1 0 0 0 0 0 0
L6_L8_rep1 1 0 0 0 0 0 0
L6_L8_rep2 1 0 0 0 0 0 0
Q5_Q7_rep1 1 0 1 0 0 0 0
Q5_Q7_rep2 1 0 1 0 0 0 0
Q3_rep1 1 0 1 0 0 0 0
Q3_rep2 1 0 1 0 0 0 0
co_40d_A 0 1 0 1 0 0 0
co_40d_B 0 1 0 1 0 0 0
co_60d_A 0 1 0 0 1 0 0
co_60d_B 0 1 0 0 1 0 0
EB_A 0 1 0 0 0 1 0
EB_B 0 1 0 0 0 1 0
H9_A 0 1 0 0 0 0 1
H9_B 0 1 0 0 0 0 1

library(edgeR)
data_filter<- count table
edgeR.dgelist = DGEList(data_filter)
edgeR.dgelist_normal = calcNormFactors(edgeR.dgelist)
CommonDisp <- estimateGLMCommonDisp(edgeR.dgelist_normal, design.matrix)
Error in glmFit.default(y, design = design, dispersion = dispersion, offset = offset, : Design matrix not of full rank. The following coefficients not estimable: methodH9

I like to know whether my design matrix is correct?
Also, I like to compare Q method against other methods, would you please help me in making contrast?

Answer: Your **method** and **batch** are confounded ( b1 -> L, Q ) and ( b2 -> co_40d, co_60d_A, EB, H9 ), same for **sample** and **batch**. Also, for b2, **method** and **sample** are the same, and it is not possible to estimate effects for redundant variables. Are **samples** technical replicates? You could drop **samples** and keep only **method**, but in any case **batch** is still confounded - meaning you can't independently estimate **batch** and **methods** effects.

Try searching for `Design matrix not of full rank`, there will be plenty of very good posts explaining the causes and how to solve it. For example:

https://support.bioconductor.org/p/68092/

https://www.biostars.org/p/199519/

https://support.bioconductor.org/p/80408/
Question: I'm currently trying to calculate the number of single-cell ATAC fragments that lie at peaks defined by a bed file.

I have a tsv file containing fragments with the following columns: chromosome, start, end, cell barcode. This file is around 10gb in size (around 250 million rows).

I also have another tsv file containing peaks with the following columns: chromosome, start, end.

What I want to do is:
1. Assign fragment to each peak on bed file
2. Output a data frame containing cell barcode, region where the fragments are assigned (chr, start, end), and also the number of fragments assigned to this region.

Here is what I tried, however R ran out of memory (I am running a desktop with 64gb memory and an i7 9700k). It also seems very slow. Maybe R is not suited for this kind of task, but currently I am most familiar with R.

In summary what I did was:
loop over fragments and filter the bed file if the fragment lies in between the region of the bed file, and then append another column to the fragments so that I have the region identifier, and finally count the unique identifier for each cell.

fragments$region <- "none"

for(i in 1:nrow(fragments)){

chr_fragment <- fragments[i,] %>% pull(chr)
start_fragment <- fragments[i,] %>% pull(start)
end_fragment <- fragments[i,] %>% pull(end)

# region containing start of fragment
extracted_region <- bed_file %>% filter(chr == chr_fragment) %>% filter(start <= start_fragment & start_fragment <=
end) %>% paste(collapse = "_")

extracted_region <- paste0("chr", extracted_region)

# append to initial fragments tsv file
fragments[i,5] <- extracted_region

paste0("finished assigning ", i, "out of ", nrow(fragments), " fragments") %>% print()
}

fragments <- rename(fragments, region = 5)

# remove fragments not in peak
fragments <- fragments %>% filter(region != "none")

# count number of fragments in region per cell
fragments <- fragments %>% mutate(dplyr::count(fragments, barcode, region))

Is it possible to manipulate this in a way that it won't take too much memory? Or is it impossible to achieve using R? Or if there is any tool that can do this (I tried to look but did not find any that suits this). It also seems to take awfully long time to perform this loop, which I don't even think will finish in weeks. Any help is greatly appreciated!

Answer: I am not sure if you have looked up for any tools for scATAC analysis. [Signac][1], part of Seurat framework is developed for chromatin data set analysis, especially scATAC.

The [FeatureMatrix][2] function does exactly what you are looking for.

FeatureMatrix(
fragments,
features,
cells = NULL,
chunk = 50,
sep = c("-", "-"),
verbose = TRUE
)

`fragments` are your fragment file and `features` are your peaks as a GRanges object.

[1]: https://satijalab.org/signac/index.html
[2]: https://satijalab.org/signac/reference/FeatureMatrix.html
Question: I have to plot some MLPA data using circos and the idea is to focus on cytobands. Circos uses the cytobands as described by UCSC/Ensembl. Looking at some genes in GeneCards I found different cytoband locations for the same gene (e. g. SCN5A). Now, I want to compare Ensembl and Entrez Gene cytoband systems and see how the plots behave. I was unable to find a file/page/table similar to UCSC goldenPath cytoBand at NCBI. Could someone point the direction, please?

Answer: Do these files from the NCBI MapView FTP site help?

wget ftp://ftp.ncbi.nih.gov/genomes/MapView/Homo_sapiens/objects/BUILD.37.3/initial_release/ideogram_9606_GCF_000001305.13_400_V1
wget ftp://ftp.ncbi.nih.gov/genomes/MapView/Homo_sapiens/objects/BUILD.37.3/initial_release/ideogram_9606_GCF_000001305.13_550_V1
wget ftp://ftp.ncbi.nih.gov/genomes/MapView/Homo_sapiens/objects/BUILD.37.3/initial_release/ideogram_9606_GCF_000001305.13_850_V1

# examine contents
# head -5 ideogram_9606_GCF_000001305.13_400_V1

1 p 36.3 0 451 1 7200000 gneg
1 p 36.2 451 682 7200000 16200000 gpos 100
1 p 36.1 682 1259 16200000 28000000 gneg
1 p 35 1259 1583 28000000 34600000 gpos 100
1 p 34.3 1583 1779 34600000 40100000 gneg
Question: I'm a Computer Science undergraduate writing my thesis which is about analyzing open source bioinformatics projects and extracting standalone modules. I need those projects to be written in Java and preferably be stable or mature. So far I've found a bunch of projects via SourceForge and I was wondering if there are any notable ones that are not hosted there. A different repository that has a bioinformatics category would be helpful too.

Here's a list of the most promising projects I have so far:

- Jenetics
- JGap
- Jmol
- Juicebox
- Picard
- RDP Classifier
- Shim
- Tassel
- The Chemistry Development Kit
- Toxfree
- TreeView
- VarScan

Thanks in advance.

Answer: - [htsjdk][1]
- [gatk][2]
- [igv][3]
- [haploview][4]

...
**Edit** and [using the seqanswer wiki][5].

[1]: https://github.com/samtools/htsjdk/
[2]: https://www.broadinstitute.org/gatk/
[3]: https://www.broadinstitute.org/igv/
[4]: https://www.broadinstitute.org/scientific-community/science/programs/medical-and-population-genetics/haploview/haploview
[5]: http://seqanswers.com/w/index.php?title=Special%3AAsk&q=[[Language%3A%3AJava]]&po=&sort_num=&order_num=ASC&eq=yes&p[format]=broadtable&p[limit]=&p[sort]=&p[offset]=&p[headers]=show&p[mainlabel]=&p[link]=all&p[searchlabel]=&p[intro]=&p[outro]=&p[default]=&p[class]=sortable+wikitable+smwtable&eq=yes
Question: Hello! everyone:

I'm new to linux, here I got a problem:

I have a file file1 like:

3
6
7
9
12

and file2 which is tab-delimited:

chr1 3052600 3052800 1 E3
chr1 3052800 3053000 2 E3
chr1 3059400 3059600 3 E3
chr1 3059600 3059800 4 E3
chr1 3059800 3060000 5 E3
chr1 3062600 3062800 6 E3
chr1 3101000 3101200 7 E3
chr1 3105000 3105200 8 E3
chr1 3105200 3105400 9 E3
chr1 3116800 3117000 10 E2
chr1 3117000 3117200 11 E2
chr1 3164800 3165000 12 E2

and I want to extract the lines in file2 which its 4-th column equal the number in file1 like below:

chr1 3059400 3059600 3 E3
chr1 3062600 3062800 6 E3
chr1 3101000 3101200 7 E3
chr1 3105200 3105400 9 E3
chr1 3164800 3165000 12 E2

I have spent several hours including wrote a very slow python script, and I searched for the oneline solution, but I got nothing!

awk -v FS="\t" 'NR==FNR{rows[$1]++;next}(substr($NF,1,length($NF)-1) in rows)' fiel1 file2

Thanks a lot for some suggestions!

Answer: simple with https://github.com/shenwei356/csvtk

csvtk grep -H -t -f 4 -P file1 file2 > result
Question: Is there a tool like bedtools shuffle which I can use to randomly shuffle a bam file? or will I have to convert my bam into a bed and then shuffle it? thanks

Answer: use https://github.com/lindenb/jvarkit/wiki/Biostar145820 to shuffle the reads with option `-n -1`
Question: Hi everyone

Is there a fast way to do this filter?

I have a huge Fasta file (sequences are short reads coming from an Illumina instrument). I have also a list of nucleotide sequences (not Fasta, just the sequences) and I want to remove from the big Fasta file all entries identical to those in the list.

My idea was simply to go down through the Fasta file and then, for every read, check all the sequences of the list. If the read matches one of the sequences then do nothing, otherwise print the read into a new file. I made this with perl but it takes ages!

The list is made up of nucleotide sequences, not IDs. It's something like this:

AACGACTACTTATCGATC

TCGGCGATATACGTAC

CCAGTTTCGGGGCTAT
....

Thanks!

Answer: linearize your fastq file into 3 columns: sequence, ID, qualities (using 'awk', or 'perl')

sort it using the sequence as the key ( unix 'sort')

sort your 2nd file of sequences ( unix 'sort')

use unix 'join' to filter out the first file and re-transform it to FASTQ using awk or perl
Question: I have a data table that I would like to run expression analysis on it but the data set is present in log2. I would like to know how can I get Raw value instead of log2?

Data: in Log2

Gene_names G_rep1.Log2 G_rep2.Log2 G_rep3.Log2
ACTG1 38.49 38.53 33.39
TUBB4B 37.61 37.31 37.35
TUBA1B 37.53 37.30 36.99
RPLP2 33.20 32.90 33.78
ACTC1 35.19 35.35 29.73

Answer: 2^x
Question: Hi all

I have a large data frame and I am trying to duplicate each column right beside the original one.

For example,

Input

V1 V2 V3
A T G

output

V1 V1 V2 V2 V3 V3
A A T T G G

Please let me know if there is a way to do this. Would it be possible to use rep() for this. Thanks in advance!

Answer: I don't know. Is this what you need?

df <- data.frame(
a = c(1,2,3,4,5),
b = c('a','b','c','d','e'),
c = c(1,'b',3,'d',5),
d = c('a',2,'c',4,'e'))
df
a b c d
1 1 a 1 a
2 2 b b 2
3 3 c 3 c
4 4 d d 4
5 5 e 5 e

via `apply()`
---------------

do.call(cbind, apply(df, 2, function(x) data.frame(x,x)))
a.x a.x.1 b.x b.x.1 c.x c.x.1 d.x d.x.1
1 1 1 a a 1 1 a a
2 2 2 b b b b 2 2
3 3 3 c c 3 3 c c
4 4 4 d d d d 4 4
5 5 5 e e 5 5 e e

via `lapply()`
-----------------

do.call(cbind, lapply(df, function(x) data.frame(x,x)))
a.x a.x.1 b.x b.x.1 c.x c.x.1 d.x d.x.1
1 1 1 a a 1 1 a a
2 2 2 b b b b 2 2
3 3 3 c c 3 3 c c
4 4 4 d d d d 4 4
5 5 5 e e 5 5 e e

Now going out for a jog. Catch you later

Kevin
Question: I am interested in adding multiple row annotations to my heatmap using the pheatmap function in R (example: https://stackoverflow.com/questions/41628450/r-pheatmap-change-annotation-colors-and-prevent-graphics-window-from-popping-up). The examples I see online all give only one annotation. What is the easy way to extend it to multiple row annotations?

Answer: from the pheatmap help:

```
# make test matrix
test = matrix(rnorm(200), 20, 10)
test[1:10, seq(1, 10, 2)] = test[1:10, seq(1, 10, 2)] + 3
test[11:20, seq(2, 10, 2)] = test[11:20, seq(2, 10, 2)] + 2
test[15:20, seq(2, 10, 2)] = test[15:20, seq(2, 10, 2)] + 4
colnames(test) = paste("Test", 1:10, sep = "")
rownames(test) = paste("Gene", 1:20, sep = "")

# define the annotation
annotation_row = data.frame(
GeneClass = factor(rep(c("Path1", "Path2", "Path3"), c(10, 4, 6))),
AdditionalAnnotation = c(rep("random1", 10), rep("random2", 10))
)
rownames(annotation_row) = paste("Gene", 1:20, sep = "")

pheatmap(test, annotation_row = annotation_row)
```
Question: Hi All.
I am using StringTie to assemble transcriptome from my RNA-Seq data. The question is that if I use refSeq as reference annotation which was download from UCSC genome browser website, does the CDS, start codon and stop codon segments in that gtf file will affect transcriptome assembly? Like would StringTie consider CDS/start codon/stop codon as a new exon, but actually these features are just parts of exons?

The gtf file looks like this:

chr1	hg19_refGene	start_codon	67000042	67000044	0.000000	+	.	gene_id "NM_032291"; transcript_id "NM_032291";

chr1	hg19_refGene	CDS	67000042	67000051	0.000000	+	0	gene_id "NM_032291"; transcript_id "NM_032291";

chr1	hg19_refGene	exon	66999639	67000051	0.000000	+	.	gene_id "NM_032291"; transcript_id "NM_032291";

chr1	hg19_refGene	CDS	67091530	67091593	0.000000	+	2	gene_id "NM_032291"; transcript_id "NM_032291";

chr1	hg19_refGene	exon	67091530	67091593	0.000000	+	.	gene_id "NM_032291"; transcript_id "NM_032291";

Answer: hi,
Plz. do not use this GTF file from UCSC. A GTF file not only has the information of individual exons (of a transcript isoform) but also of different transcripts (that originate from a particular gene).
You would notice that the `gene_id` and `transcript_id` are same in the above GTF file. So any transcript assembler you use (like StringTie) would not be able to infer the transcript <-> gene relationship.

See this GTF structure from Ensembl -

1	protein_coding	exon	874655	874840	.	+	.	gene_id "ENSG00000187634"; transcript_id "ENST00000455979"; exon_number "1"; gene_name "SAMD11"; gene_biotype "protein_coding"; transcript_name "SAMD11-004"; exon_id "ENSE00002715021";
1	protein_coding	CDS	874655	874840	.	+	2	gene_id "ENSG00000187634"; transcript_id "ENST00000455979"; exon_number "1"; gene_name "SAMD11"; gene_biotype "protein_coding"; transcript_name "SAMD11-004"; protein_id "ENSP00000412228";
1	protein_coding	exon	876524	876686	.	+	.	gene_id "ENSG00000187634"; transcript_id "ENST00000455979"; exon_number "2"; gene_name "SAMD11"; gene_biotype "protein_coding"; transcript_name "SAMD11-004"; exon_id "ENSE00003477353";
1	protein_coding	CDS	876524	876686	.	+	2	gene_id "ENSG00000187634"; transcript_id "ENST00000455979"; exon_number "2"; gene_name "SAMD11"; gene_biotype "protein_coding"; transcript_name "SAMD11-004"; protein_id "ENSP00000412228";

Hope this is clear. Plz use GTF from [Ensembl][1] or [Gencode][2]

[1]: ftp://ftp.ensembl.org/pub/release-84/gtf/homo_sapiens/
[2]: http://www.gencodegenes.org/releases/current.html
Question: Hello All,

I was curious if anyone is aware of any tools or approach that can convert hard-masked genome to soft-masked format. For example,

genome: *....ATGCATGCATGC......*

conversion required: *....ATGCNNNNATGC......* to *....ATGCatgcATGC......*

Regards,
B

Answer: I would first get the coordinates of Ns in the hard-masked genome and output them as a bed file. Here's an example using seqkit.

```
seqkit locate --bed -rPp "N+" hardmasked.fasta > N_coords.bed
```

You can then use bedtools to soft mask the non-masked genome.

```
bedtools maskfasta -soft -fi unmasked.fasta -bed N_coords.bed -fo softmasked.fasta
```
Question: Hi Biostars,

I got two output files (*genes.results and *isoforms.results) from aligning reads to the Trinity denovo assembly(Trinity.fasta) with Bowtie2 and expected read counts with RSEM. Out of two output files, which one should I feed to DESeq2 for identifying DEGs.

Any suggestions will be appreciated. Thanks in advance.

Answer: That means if I want gene level LFC I need to input *genes.results and for transcript level LFC *isoforms.results file should be used. Thank you.
Question: Hi,

Assume a table as below:

X =

col1 col2 col3
row1 "A" "0" "1"
row2 "B" "2" "NA"
row3 "C" "1" "2"

I select combinations of two rows, using the code below:

pair <- apply(X, 2, combn, m=2)

This returns a matrix of the form:

pair =

[,1] [,2] [,3]
[1,] "A" "0" "1"
[2,] "B" "2" NA
[3,] "A" "0" "1"
[4,] "C" "1" "2"
[5,] "B" "2" NA
[6,] "C" "1" "2"

I wish to iterate over pair, taking two rows at a time, i.e. first isolate `[1,]` and `[2,]`, then `[3,]` and `[4,]` and finally, `[5,]` and `[6,]`. These rows will then be passed as arguments to regression models, i.e. `lm(Y ~ row[i]*row[j])`.

I am dealing with a large dataset. Can anybody advise how to iterate over a matrix two rows at a time, assign those rows to variables and pass as arguments to a function?

Thanks,
S ;-)

**Edit:** In response to the comments, I should specify that my problem concerns SNP and expression data where I aim to do a pairwise multiple regression analysis (first order regression) in order to assess any possible SNP-SNP interactions that may effect the expression phenotype.

Answer: A simple idiom for "all pairs of a matrix" M of dim(N, C)

for( i in 1:(N-1) ){
for( j in (i+1):N ){
print(paste(i,j))
}
}

I've had to do this eQTL analysis. For any non-trivial number of genotypes or expression measurements, the naive approach (e.g. building a new lm for each comparison) is too slow. If PHENO is a matrix (not data.frame) of expression values and GENO the set of genotypes for a single locus at all of the samples in PHENO, you want something like

for each GENO
model = lm( t(PHENO)~GENO )

Extend the linear model to match whatever interactions you want to test. You can extract the p-values of coefficients of a solved linear model with

calculate.lm.pval=function(linear.model){
# Given a solved linear model, return the p-values of the coefficients
n.obs = dim(linear.model$residuals)[1]
n.models = dim(linear.model$coefficients)[2]
cc = matrix(coef(linear.model), nrow=linear.model$rank, ncol=n.models)
se = calculate.linear.model.se(linear.model)
t.stat = cc/se
rdf = n.obs - linear.model$rank
dm=data.frame( 2*pt(abs(t.stat), rdf, lower.tail = FALSE) )
names(dm) = as.vector(dimnames(linear.model$coefficients)[[2]])
rownames(dm) = dimnames(linear.model$coefficients)[[1]]
dm
}
Question: Dear Users,

I would like to know , that given I have only one tumor and its normal type and I have the sets of germline and somatic mutations with the frequency listed by VarScan2, how can I use this information along with the region information of this mutation corresponding to a gene , to outline the clonal and sub-clonal populations of mutation in the tumor. Is there any method that can help me generate a model which can help me understand which mutations are occupied in the entire tumor population and which less. This would help me understand all the sub clones of the tumor. In other way it will help me understand how the mutations are categorized in the entire tumor mass and also inform me to what extent this mutation is having a stake in the tumor. This classification helps to reconstruct the tumor fate and its evolution and will also enable me to list out the potential driver mutations and passenger mutations concerned with that tumor. Is there any tool that can help me do this? Most of the tools work on multiple samples. I would like to have some suggestions on these lines.

Answer: We built the sciClone package for exactly this purpose: https://github.com/genome/sciclone

It takes inputs of somatic mutations, with readcounts and VAFs, and uses that information to infer subclonal populations in heterogeneous tumors. It also gives you some nice visualization options.