cannin/biostars_qa huggingface dataset:
Question: Why is it so difficult to make things in ggplot2 , i like the way it helps in customisation but the curve is steep nevertheless

Here is my sample dataframe

df <- gene HSC CMP
ENSG00000158292.6 1.8102636 2.456869
ENSG00000162496.6 2.6796705 6.203838
ENSG00000117115.10 3.4509115 5.555739
ENSG00000159423.14 3.6809277 5.063446
ENSG00000053372.4 5.7089974 6.851090

If i have plot a boxplot i can simply write this `boxplot(df[,-1],col=c("red","blue"))`

I get a boxplot but when im trying with ggplot2 im having difficult time

ex <- melt(df, id.vars=c("HSC", "CMP"))
ggplot(data = ex,
aes(x = CMP, y = HSC)) +
geom_boxplot()

I get a single boxplot what i want is i get a box plot for HSC and CMP which i got when i use simple base R boxplot .

Any help or suggestion would be highly appreciated with my ggplot2 code

Answer: options(stringsAsFactors = F)
df= read.csv("test.txt", sep="\t")
library(reshape2)
library(ggplot2)
df_melt=melt(df,id.vars="gene")
ggplot(df_melt, aes(variable,value)) +
stat_boxplot(geom="errorbar", width=.5)+
geom_boxplot(aes(fill=variable))+
theme_bw()+
theme(axis.title.x=element_blank(), axis.title.y=element_blank())+
stat_summary(fun.y=median, colour="red", geom="line", aes(group = 1))+
geom_jitter(position = position_jitter(0.2))

Input:

> df
gene HSC CMP
1 ENSG00000158292.6 1.810264 2.456869
2 ENSG00000162496.6 2.679670 6.203838
3 ENSG00000117115.10 3.450912 5.555739
4 ENSG00000159423.14 3.680928 5.063446
5 ENSG00000053372.4 5.708997 6.851090
Question: Hi there,

This is my first time performing chip-seq analysis. I would like to use SICER to identify enriched domains from histone modification data. One of the input parameters that SICER requires is the fragment size. From the manual:

> Fragment size: is for determination of the amount of shift from the beginning of a read to the center of the DNA fragment represented by the read. FRAGMENT_SIZE=150 means the shift is 75.

How can I calculate the fragment size?

Answer: There is a script inside SICER download folder, called fragment-size-estimation.sh , you can try to use it to estimate the fragment size. I would recommend you to read SICER users google group. There are some interesting posts like: https://groups.google.com/forum/#!searchin/sicer-users/fragment/sicer-users/pWF0GmYYAMM/IISEiBArg18J talking about fragment size estimation issue.

Hope it helps.
Question: Hi,

I am trying to work with "minet" package in R. I have to import a large table of expression values in R. in minet tutorial it has said that data must be entered as data frame using this command: `data(syn.data)`. I do not know in which format I have to import my data or what file extension should I use.

Can anybody suggest me any solution?

Regards
Nazanin

Answer: Hi

My problem was solved. The input I used was a compendium of gene expression data, i.e. the rows were gene names and the columns correspond to samples. Before using minet I had to transpose this matrix in excel, because minet needs genes as columns and samples as rows. After reading this table via read.table command, minet works properly.

Let me know if you have any question

Good luck

Nazanin
Question: I have a data frame in which the first column contains unique identifiers in the form of Arabidopsis AGI transcript numbers.

target_id	est_counts	tpm
AT1G01010.1	236	26.3
AT1G01020.2	55.0	12.2
AT1G01020.6	0	0
AT1G01020.1	25.2	4.8
AT1G01020.4	5.8	0.2
AT1G01020.5	45.5	8.8
AT1G01020.3	3.5	0.5
AT1G01030.2	13.25	1.3
AT1G01030.1	17.75	1.7

This table has data for three genes; AT1G01010, AT1G01020 and AT1G01030. However, two genes, AT1G01020 and AT1G01030 have multiple transcripts as indicated by the number to the right of the decimal. I would like to collapse the above table into three entries, one entry for each gene in which each column will contain the total for all transcripts of that gene. So, the resulting data frame will look like this:

target_id	est_counts	tpm
AT1G01010	236	26.3
AT1G01020	135.0	26.5
AT1G01030	31	3.0
Here, the values for est_counts and tpm for all transcripts of a single gene have been summed. I have seen things like lapply(.sd), but I don't think that will work because the AGI numbers are not exactly equal and my real table has thousands of different AGI numbers.

Answer: Clean up `target_id` column, then aggregate other multiple columns at once, see:

# reproducible data
df1 <- read.table(text = "target_id est_counts tpm
AT1G01010.1 236 26.3
AT1G01020.2 55.0 12.2
AT1G01020.6 0 0
AT1G01020.1 25.2 4.8
AT1G01020.4 5.8 0.2
AT1G01020.5 45.5 8.8
AT1G01020.3 3.5 0.5
AT1G01030.2 13.25 1.3
AT1G01030.1 17.75 1.7", header = TRUE, stringsAsFactors = FALSE)

# we can imagine they are filenames and drop "extensions":
df1$target_id_clean <- tools::file_path_sans_ext(df1$target_id)

# then aggregate and sum on multiple columns
aggregate(.~target_id_clean, df1[-1], sum)
# target_id_clean est_counts tpm
# 1 AT1G01010 236 26.3
# 2 AT1G01020 135 26.5
# 3 AT1G01030 31 3.0
Question: This is somewhat of a continuation to my previous question [here][1]. Apparently, the problem with the paired-end artificial reads is that they have the following orientation:

```
-----> read 1 -----> read 2
--------------------------------- DNA fragment
```

Is there a Bowtie 2 configuration that will allow for these reads to be aligned as concordant pairs as opposed to discordant ones?

[1]: /p/104086/#104104

Answer: From the manual:

> The upstream/downstream mate orientations for a valid paired-end alignment against the forward reference strand. E.g., if `--fr` is specified and there is a candidate paired-end alignment where mate 1 appears upstream of the reverse complement of mate 2 and the fragment length constraints (`-I` (The minimum fragment length for valid paired-end alignments) and `-X` (The maximum fragment length for valid paired-end alignments)) are met, that alignment is valid. Also, if mate 2 appears upstream of the reverse complement of mate 1 and all other constraints are met, that too is valid. `--rf` likewise requires that an upstream mate1 be reverse-complemented and a downstream mate2 be forward-oriented. `--ff` requires both an upstream mate 1 and a downstream mate 2 to be forward-oriented. Default: `--fr` (appropriate for Illumina's Paired-end Sequencing Assay).

So.. `--ff`?
Question: Are there some faster alternatives to BLAST (specifically nucleotide BLAST)? I like that I can search through all GenBank+EMBL+DDBJ+PDB+RefSeq sequences (nt collection), but I feel like there must be a faster way. If I wanted to identify thousands or millions of sequences, it's somewhat inefficient.

Answer: I guess what I was really asking for is a metagenomic classifier. There are a few of those out there:

- Kraken: https://ccb.jhu.edu/software/kraken/
- GOTTCHA: http://lanl-bioinformatics.github.io/GOTTCHA/
- CLARK: http://clark.cs.ucr.edu/
- MetaPhlAn: http://huttenhower.sph.harvard.edu/metaphlan

Sure, these are not exactly same as BLAST, but if you need to quickly classify a lot of reads, these tools will do that.
Question: Hi all,

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

Thanks,
Deepthi

Answer: Another possibility, from BBMap:

readlength.sh in=reads.fq out=histogram.txt

The default is 10bp bins with a max of 80kbp, but those can be configured (run the shellscript with no arguments for details). It's very fast, and handles fastq/fasta/sam/bam; raw/gzip/bzip2.
Question: Good Morning, I'm sorry to disturb, but as anyone done these volcano plot? I'm only able to do the traditional one, i'm kind knew too these field.
Thanks Ana Sofia Moreira

Answer: **Edit 28th March, 2019:**

this is now a Bioconductor package: https://www.biostars.org/p/335751/

--------------------
-------------

That looks like an up-side down volcano plot, which is unusual. Why not just do a regular one?

You can use my EnhancedVolcano function on my GitHub page (in part developed with a Junior Bioinformatician): https://github.com/kevinblighe (https://github.com/kevinblighe/EnhancedVolcano)

Apart from anything else, it makes the use of `ggrepel` in order to fit as many gene labels as possible.
Question: Dears,

I want to know whether there is difference between different genome fasta files like `Mus_musculus.GRCm38.dna.toplevel.fa` file downloaded from ENSEMBL and `chrom_mm10.fa` file downloaded UCSC. Whether we can annotate the alignment bam files (using UCSC genome as reference) with ENSEML annotations.

Any comments and explanations are appreciated.

Answer: I would suggest reading the README file for the data. [Here is the Ensembl one.](ftp://ftp.ensembl.org/pub/current_fasta/mus_musculus/dna/README)

Generally, it is recommended that you use reference annotation and genome files from the same source.
Question: Following up on my [previous question](https://www.biostars.org/p/209854/), I have inspected the length distribution of the first exons of spliced transcripts in GENCODE26.

I still find it surprising that some are tiny, often starting at the translation start site, while the transcript is not annotated as "processed". A list of all first spliced exons shorter than 4 nt can be found at the bottom of the following page and the source code generating it is also available [on GitHub](https://github.com/charles-plessy/gencode-feature-lengths).

Is there a better explanation than just imperfections in processing pipelines or annotations ?

The reason I ask is that I would like to know if there is a lower limit of first exons at the time of RNA splicing. (And the reason behind the reason is that I have a sequence library from which I can infer the length of first exons, and I would like to test if the results are reliable).

![Length of first spliced exons](https://raw.githubusercontent.com/charles-plessy/gencode-feature-lengths/960c0b25eae5c52bf07dc8c6c3ea3862174b0aa1/gencode-feature-lengths_files/figure-html/first-microexons-1.png)

Answer: This is an interesting analysis of the GENCODE v26 geneset, that highlights potential idiosyncrasies in the annotation process. There are however some important points to consider when interpreting these results.

While the current geneset contains first 'exons' with less than 20nt, these do not necessarily represent the full-length of the exon. Indeed, these first exons may not actually represent the transcriptional start of the associated transcript. This is because manually annotated GENCODE transcripts are extended to the limit of the supporting evidence, which typically consists of ESTs and cDNAs. If this evidence doesn't extend to the actual TSS, then the annotated transcript will essentially be truncated at the 5' end. This appears to be the case for most of the 'micro first exons' identified here, as they coincide with longer exons in other transcript isoforms for the same gene.

Although the annotation of such short exons may be accurate, their actual value within the geneset may be limited. It may be best therefore to exclude such transcripts from any analysis of RNA splicing. On behalf of GENCODE, we welcome any comments or questions, which can be sent to either [HAVANA][1] or [GENCODE][2].

[1]: mailto:havana-help@sanger.ac.uk
[2]: mailto:gencode-help@sanger.ac.uk
Question: I would like to find the count of homozygous major(0|0),homozygous minor(1|1) and heterozygous(1|0,0|1) alleles for each position of a chromosome.I would like to form a table like this:

**POS 00 1 0/0 1 11**

Pos1 n1 n2 n3

Pos2 n4 n5 n6

where n1...n6 are the counts .

See the example file here [1000Genomes][1]

[1]: ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz

Any suggestions to try to get this information would be much appreciated

Thanks in Advance

Answer: An `awk` solution:

$ zcat input.vcf.gz|awk -v OFS="\t" '$0 !~ "^#" {hom_ref = 0; hom_alt = 0; het = 0; for(i=10;i<=NF;i++) { if($i ~ /0\|0/) hom_ref++; else if($i ~ /1\|1/) hom_alt++; else het++; } print $1, $2, hom_ref, hom_alt, het}'

fin swimmer
Question: I'm still a beginner and using [Illumina Platinum Genomes as a toy dataset][1]. Google supplies the files in VCF, but I need them annotated with rsIDs.

I've tried GATK VariantAnnotator and bcftools annotate [as suggested by this post I found via Google][2] and [this Biostars post][3], but I'm not having any luck.

When the command line finishes running, the output doesn't have the rsIDs filled in and no obvious errors pop up. I'm kinda stumped.

Any help is appreciated.

[1]: https://cloud.google.com/genomics/docs/public-datasets/illumina-platinum-genomes
[2]: https://gist.github.com/obenshaindw/99440ac43de07548e453
[3]: https://www.biostars.org/p/160021/

Answer: Hey it's quite some time ago but if anyone else is having a problem I just wanted to say following command worked for me:

bcftools annotate -a /data/references/hg19/pipe/dbsnp138/00-All.vcf.gz -c ID -o samtools_annotated.vcf.gz samtools.vcf.gz

The thing to look out for is, I think it doesn't work if your reference has numeric chromosome column while your VCF file has it as with 'chr' added to is, *e.g.* "chr1".
Question: I have a fastq file. I want to have k-mers and their abundance numbers.
can you give me some tools to do this!
thank you

Answer: If you're working with small values of k, I wrote a command-line k-mer counter called `kmer-counter` that will output results in a form that a Python script can consume:

https://github.com/alexpreynolds/kmer-counter

You can grab, build and install it like so:

$ git clone https://github.com/alexpreynolds/kmer-counter.git
$ cd kmer-counter
$ make
$ cp kmer-counter /usr/local/bin

Once the binary is in your path, you might use it in Python like so:

k = 6
fastaFile = '/path/to/some/seqs.fa'
kmerCmd = 'kmer-counter --fasta --k=%d %s' % (k, fastaFile)
try:
output = subprocess.check_output(kmerCmd, shell=True)
result = {}
for line in output.splitlines():
(header, counts) = line.strip().split('\t')
header = header[1:]
kmers = dict((k,int(v)) for (k,v) in [d.split(':') for d in counts.split(' ')])
result[header] = kmers
sys.stdout.write("%s" % (str(result)))
except subprocess.CalledProcessError as error:
sys.stderr.write("%s" % (str(error)))

Given example FASTA like this:

>foo
TTAACG
>bar
GTGGAAGTTCTTAGGGCATGGCAAAGAGTCAGAATTTGAC

For k=6, you would get an iterable Python dictionary like this:

{'foo': {'TTAACG': 1, 'CGTTAA': 1}, 'bar': {'GTTCTT': 1, 'AGAACT': 1, 'GAGTCA': 1, 'ATGGCA': 1, 'GAACTT': 1, 'ATTCTG': 1, 'CTAAGA': 1, 'CTTCCA': 1, 'ATTTGA': 1, 'GGAAGT': 1, 'AGGGCA': 1, 'CCTAAG': 1, 'CTCTTT': 1, 'AATTTG': 1, 'TCTGAC': 1, 'TTTGCC': 1, 'CTTAGG': 1, 'TTTGAC': 1, 'GAAGTT': 1, 'CCCTAA': 1, 'AGAATT': 1, 'AGTCAG': 1, 'CTGACT': 1, 'TCTTAG': 1, 'CGTTAA': 1, 'GTGGAA': 1, 'TGCCAT': 1, 'ACTCTT': 1, 'GGGCAT': 1, 'TTAGGG': 1, 'CTTTGC': 1, 'TGGAAG': 1, 'GACTCT': 1, 'CATGCC': 1, 'GCAAAG': 1, 'AAATTC': 1, 'GTCAAA': 1, 'TGACTC': 1, 'TAGGGC': 1, 'AAGTTC': 1, 'ATGCCC': 1, 'TCAAAT': 1, 'CAAAGA': 1, 'AACTTC': 1, 'GTCAGA': 1, 'CAAATT': 1, 'TAAGAA': 1, 'CATGGC': 1, 'AAGAAC': 1, 'AAGAGT': 1, 'TCTTTG': 1, 'TTCCAC': 1, 'TGGCAA': 1, 'GGCAAA': 1, 'AGTTCT': 1, 'AGAGTC': 1, 'TCAGAA': 1, 'GAATTT': 1, 'AAAGAG': 1, 'TGCCCT': 1, 'CCATGC': 1, 'GGCATG': 1, 'TTGCCA': 1, 'CAGAAT': 1, 'AATTCT': 1, 'GCATGG': 1, 'ACTTCC': 1, 'TTCTTA': 1, 'GCCATG': 1, 'GCCCTA': 1, 'TTCTGA': 1}}

You can use standard Python calls to manipulate this dictionary object and get sums of counts per record, for sequence, etc. which seems to answer your question. Please feel free to clarify what you're looking for if this object representation is not clear.

For a fully-fleshed out demonstration, see: https://github.com/alexpreynolds/kmer-counter/blob/master/test/kmer-test.py

If you are working with large values of k, then you might look at other tools, like Jellyfish. They use data structures to give an approximate count for large values of k that are not tenable for strict word counts.
Question: hi,

I have a list of miRNAs IDs from different organisms how I can take sequence of mature miRNAs for all my IDs?

for example

lgi-mir-34

mdo-mir-193

mdo-mir-21

mdo-mir-302d

mml-mir-18

Answer: From miRBase [using this file][1].

[1]: ftp://mirbase.org/pub/mirbase/CURRENT/mature.fa.gz
Question: I am running IQ-TREE on a dataset of 160 samples with on average 2.1 millions reads on a genome of 310 Mb.

However, the software has been computing for one week now, and progress in the standard output are very slow.

The last segment of the standard output are:

Create initial parsimony tree by phylogenetic likelihood library (PLL)... 101.349 seconds Perform fast likelihood tree search using GTR+I+G model... Estimate model parameters (epsilon = 5.000) Perform nearest neighbor interchange... Optimizing NNI: done in 37959.5 secs using 99.58% CPU Estimate model parameters (epsilon = 1.000)
1. Initial log-likelihood: -28771040.770
2. Current log-likelihood: -28771036.989
3. Current log-likelihood: -28771033.633
4. Current log-likelihood: -28771030.718
5. Current log-likelihood: -28771028.173
6. Current log-likelihood: -28771025.940
7. Current log-likelihood: -28771023.982
8. Current log-likelihood: -28771022.265
9. Current log-likelihood: -28771020.758
10. Current log-likelihood: -28771019.437
11. Current log-likelihood: -28771018.277
12. Current log-likelihood: -28771017.259 Optimal log-likelihood: -28771016.360 Rate parameters: A-C: 0.99469 A-G: 2.56577 A-T: 1.06937 C-G: 0.76029 C-T: 2.56126 G-T: 1.00000 Base frequencies: A: 0.254 C: 0.244 G: 0.245 T: 0.256 Proportion of invariable sites:
0.000 Gamma shape alpha: 3.383 Parameters optimization took 12 rounds (9957.005 sec) Time for fast ML tree search: 101193.759 seconds

NOTE: ModelFinder requires 214666 MB RAM! ModelFinder will test up to 484 DNA models (sample size: 4006362) ... No. Model -LnL df AIC AICc BIC 1 GTR+F 28847670.292 339
57696018.584 57696018.642 57700494.535 2 GTR+F+I 28847652.995 340 57695985.990 57695986.047 57700475.144 3 GTR+F+G4
28771009.772 340 57542699.545 57542699.602 57547188.699 4 GTR+F+I+G4 28771010.076 341 57542702.152 57542702.210 57547204.509 5 GTR+F+R2 28587738.467 341 57176158.933 57176158.992
57180661.291

Would anyone know if ModelFinder will effectively test 484 DNA models?

It took days to output these five tests. Is there any shortcut for the estimation of the right substitution model?

Answer: > I am running IQ-TREE on a dataset of 160 samples with on average 2.1 millions reads on a genome of 310 Mb.

Not quite sure what you mean based on the above explanation, but log-likelihood values tell me this is a monstrous alignment. Nothing will be done fast with large alignments.

If on top of that you have insufficient RAM and the program needs to swap a lot, this could take months. I may be able to offer better advice if you explain in greater detail what your exact alignment is, but either way it sounds like you need to rethink your strategy. That means starting from a smaller alignment, or using a generic substitution model with an alignment you have.
Question: I am running BEAGLE 4.1 on data set with three groups. I ran a .vcf that has all three groups' data combined, which went without a hitch and I have great useable data. However, I also split up the main .vcf into three .vcfs for each group using R, to impute them individually. This is where I encountered my "error" - while one group (the largest) went normally, the imputation on the second group always stops early. The command prompt, when I check it, has the same line repeated as far up as I can scroll, which reads:

at dag.MergeableDag.similar

The run log does not have this line in it, nor does it list any error; rather, it just lists windows and iterations up to the last one it ran with no sign of why it cut off. BEAGLE *does* generate an output .vcf, which is truncated somewhere early on. Every time I run BEAGLE, the exact spot it ends on differs, though the output file is always truncated with the same repeated line in the terminal. Often (but not always), the last line in the imputed .vcf is the last SNP in a chromosome. I have looked at these lines in the input .vcf and noticed nothing special about them compared to others (no missing data or strange formatting).

I have so far tried restarting the computer several times, closing all processes, updating Java, and deleting/re-downloading the Beagle .jar, but the error persists. I am running Windows 7 Professional.

Has anyone had a similar error? What does it mean?

Thank you in advance!

**Edit:** It was pointed out in an email from Dr. Browning that this was likely a stack overflow error. I had not initially been using the -Xss option, but I added `-Xss5m` to my command prompt at his suggestion and the run worked.

Answer: Adding OP's solution to close question with an accepted answer:

----
**Edit:** It was pointed out in an email from Dr. Browning that this was likely a stack overflow error. I had not initially been using the -Xss option, but I added `-Xss5m` to my command prompt at his suggestion and the run worked.

----
https://www.biostars.org/u/44123/ : If you could add this as an answer too, I'll delete my post so you're the one answering it, not me.
Question: Hi everyone,

I installed the modeller program on my computer. I'm running a Lubuntu 32-bit system. I have python installed. I can check that by running `python --version` on a terminal. However whenever I run a script that has modeller functions in it, like executing this script: `python model-default.py > model-default.log` it's giving me this result:

Traceback (most recent call last):
File "model-default.py", line 4, in
from modeller import * # Load standard Modeller classes
ImportError: No module named modeller

What should I do in order to be able to import modeller inside my python scripts? Also is there an easy API for modeller so that I can run it in java? Thank you.

Answer: When you install MODELLER it warns you to add a few folders to the PYTHONPATH and LD_LIBRARY_PATH environment variables. Locate your MODELLER installation to find out the folders you need. You can use the command: locate _modeller.so.

Our .bashrc scripts have the following in them to load the folder at login.

if [ -z $PYTHONPATH ]
then
export PYTHONPATH="/folder/science/modeller/lib/x86_64-intel8/python2.5/:/folder/science/modeller/modlib/"
else
export PYTHONPATH="${PYTHONPATH}:/folder/science/modeller/lib/x86_64-intel8/python2.5:/folder/science/modeller/modlib/"
fi

if [ -z $LD_LIBRARY_PATH ]
then
export LD_LIBRARY_PATH="/folder/science/modeller/lib/x86_64-intel8"
else
export LD_LIBRARY_PATH="${LD_LIBRARY_PATH}:/folder/science/modeller/lib/x86_64-intel8"
fi