cannin/biostars_qa huggingface dataset:
Question: Hello,

Please help me with understanding my boxplot.

Here is the boxplot, I got for my RNA-Seq data.

![Box Plot][1]

My data is

```r
head(rawCountTable)
con-1 con-2 a-1 a-2 b-1 b-2 c-1 c-2 d-1 d-2
ENSMUSG0000000000 0 0 0 0 0 0 0 0 0 0
ENSMUSG00000000028 854 937 1143 1029 912 856 809 754 513 520
ENSMUSG00000000031 822918 817451 716860 691396 763705 829274 838094 819312 717935 730879
```

The code for Boxplot is below:

```r
pseudoCount = log2(rawCountTable + 1)
df = melt(pseudoCount, variable.name = "Samples",
value.name = "count") # reshape the matrix
df = data.frame(df, Condition = substr(df$Samples, 1, 4))
```

Here is my code for the density plot.

```r
ggplot(df, aes(x = count, colour = Samples, fill = Samples)) +
ylim(c(0, 0.17)) +
geom_density(alpha = 0.2, size = 1.25) +
facet_wrap(~ Condition) + theme(legend.position = "top") +
xlab(expression(log[2](count + 1)))
```

The density Plot is

![Density Plot][2]

So, my question is I want to know how to interpret these plots? How is my data quality? If you can recommend me an article about understanding these plots and assess my data, I would appreciate it.

Thank you very much!

[1]: https://i.ibb.co/TtQRfw6/Boxplot.png
[2]: https://i.ibb.co/gMdSsqN/Density-Plot.png

Answer: I don't think those plots are necessarily too informative about quality. If you want a general idea about the quality of the sequencing reads, use a program like [FastQC](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The alignment statistics from your aligner will then give you a good idea of the complexity of your library. If you plan on running differential expression on your data, you can generate PCA and heatmap plots, which will be a good first indicator of replicate concordance, and from those plots you can sometimes start seeing the difference between conditions. The [DESeq2](https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html) is a good resource for making these plots.
Question: Does anyone know how/where to access MeDUSA for methyl-seq studies? The paper describing it is [here][1].

But the link to download it does not function (http://www2.cancer.ucl.ac.uk/medicalgenomics/medusa/).

Has MeDUSA been discontinued?

**EDIT**: MeDUSA is still offline (as of 1/22/2015).

[1]: http://www.ncbi.nlm.nih.gov/pubmed/23587164

Answer: It's finally back online [here][1]. The full reply I got (they attached the source code, which is the "this file" being referred to) from the authors is below:

> I have added this file as version 2.1 to the medusa webpage at https://www.ucl.ac.uk/cancer/medical-genomics/medusa
Supporting files are available at the Gigascience webpage (http://gigadb.org/dataset/100035) as indicated in https://www.ucl.ac.uk/cancer/medical-genomics/medusaproject
The two medical genomics pages have just been published so changes should be effective online in less than 2 hours.
Let me know if there is any problem with the webpages.

[1]: https://www.ucl.ac.uk/cancer/medical-genomics/medusa
Question: Hi Guys,

I am about to start working on a biomarker discovery project using transcriptomics and I basically have no experience in biomarker discovery (although I have experience with rna-seq data analysis).

Also from theoretical point of view, I am informed of how biomarker discovery using rna-seq works.

My problem is that I don’t know the state of the art bioinformatics tools that are being used in this field. Therefore, I looked for a tutorial or a step by step workflow that I can use to quickly catch up the current best practice for transcriptomics based biomarker Discovery but I can’t find a good one.

I really appreciate if someone can introduce a good tutorial or a step by step workflow.

Thanks in advance

Answer: I never came across such tutorial...

If you are familiar with RNAseq data analysis for differential expression, you can follow from there.

What I often see is:

- Exploratory data analysis and quality control (PCA, MDS, sanity checks etc);
- Differential gene/transcript analysis (interesting to select the most differentially expressed genes in average, not necessarily the most discriminant ones);
- Depending on your sample size and what type of dependend variable you have (categorical two level, categorical multilevel, ordinal, numeric, etc.), you can train some models, e.g, penalized regression (lasso, ridge, elastic net), random forests, support vector machines, kNN etc;
- Validate on external datasets (public) or in new samples (or test set from your own experiment, not used in any steps before).

There are more things to do, but depends on study design, data quality, biological phenomenon and sample size and study objectives.

I suggest start reading good papers from reputable journals with similar study design. Also, there are some good reviews on transcriptomic biomarkers and rnaseq biomarkers or for classification purposes.
Question: Have recently signing up to Galaxy https://usegalaxy.org/root

and used it to help me - non-bioinformatician - to do some work. At Galaxy, they provide a Biostar forum that allows users to post questions regarding data processing using Galaxy. I mistakenly thought that was the same Biostar I have been using and thus was expecting to get some quick feedback for my question. It turns out they are likely two different, un-connected forums. Does any user here at Biostar-bioinformatics explained use or know about Biostar-Galaxy explained ? I am asking this because I am not sure if I am allowed to post the same question twice, if the two forums are somehow linked. Any thought? Thank you.

Answer: The galaxy site is essentially a galaxy-specific version of this site. Biostars is an open source project that's used by a few different sites (e.g., Bioconductor) for their support/help forums.

As a general rule, cross-posting on multiple forums is discouraged, regardless of whether the forums are connected/related. For example, if you post the same question here and on SEQanswers, which is in no way connected to this site, I'll see that and flag your post as being cross-posted. Similar things would likely happen with the Bioconductor/Galaxy/Neurostar/etc. sites as well. The reason is that you end up doubling the communities efforts when you cross-post, which isn't very nice.
Question: I have been asked to recommend introductory books and resources to R and Bioconductor. My problem is just, I never read a book to learn R or Bioconductor, so I have no experience with this and cannot recommend one. I am interested in mainly introductory books, possibly targeting various groups of readers (computer scientists, molecular biologists, (bio-)statisticians), any recommendation appreciated.

For example, I used the following resources:

The R-manuals, especially the R intro

There are also a lot of contributed documents there on the R web site, but I didn't
use.

If a package from Bioconductor interests me, I read the package vignette.

I read the Bioconductor mailing list, that helps to see what other people use.

I have the "Venebles, Ripley. S Programming" book, that is hardly introductory.

Which books did you find helpful or completely useless to learn R/Bioconductor? For example: R Programming for Bioinformatics looks promising, anybody read it?

Or do you share my reluctance towards R-books and prefer online resources?

Answer: I love plotting:

http://books.google.com/books?id=F_hwtlzPXBcC&printsec=frontcover&dq=ggplot2&hl=en&sa=X&ei=D7I1T_7oGcPmiALezJCUCg&ved=0CDIQ6AEwAA#v=onepage&q=ggplot2&f=false

ggplot is an amazing [R] package. I use it every day.
Question: Hello Biostars,

Today, I'm doing network analysis using Cytoscape.
I installed 'stringApp' plugin for getting STRINGDb network.
After adding 'gene symbol' list in it, I ran strinApp. As a result, got node, edge and network text files.
The node file I got contains those genes also which I never entered in input!

Has anyone faced such problem? or is it something very obvious?

Please help me out.
Thanks in advance!

Answer: In Cytoscape, the stringApp installs into the network search bar at the top of the Control Panel. I think this is where you are pasting/typing your list of gene symbols, which is correct. Before clicking on search (the magnifying glass icon), however, you will want to click on settings (the three line icon). Here you can specify how many "additional interactors" you want. Set it to zero and then you'll just get interactions between you original input.
Question: Hi, I have a list of variants with varints id like 1_98929_A_G_b37 and repsonding GENE id like ENSG00000237094.7. Is there any R package or on-line tools that can convert these kind of hg 19 SNp id and ensenble id into rs id and ordianry gene names?
Any ideas is appreciated

Answer: If you've got the list of variants in VCF format, I'd recommend using the VEP tool
[http://www.ensembl.org/info/docs/tools/vep/index.html][1]

It enables further annotation with gene name symbols and dbgap accessions etc.

[1]: http://www.ensembl.org/info/docs/tools/vep/index.html
Question: Is there a method or package in R to generate network or perform network analysis using R

Answer: My tutorial appears to be a perfect fit:
https://www.biostars.org/p/285296/

If you already have your own data matrix, then start from *Step 2*

Kevin
Question: Dear all,

I have a multi-sample (>100) VCF file and a list of SNPs with CHR and POS for which the sample names and their genotypes need to be extracted.

For example, i have a SNP chr1:23455. I would like to know how many and which samples have homozygous alternate allele, heterozygous allele and homozygous reference allele.

Is there any existing tool that provides the above summary?

Many thanks!!!

Answer: use bioalcidaejdk http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html and the following script:

```
// all the types of Genotypes see https://samtools.github.io/htsjdk/javadoc/htsjdk/htsjdk/variant/variantcontext/GenotypeType.html
final List all_types = Arrays.asList(GenotypeType.values());
//print the header
println("CHROM\tPOS\tREF\tALT\t"+all_types.stream().map(T->"count."+T.name()+"\tsample."+T.name()).collect(Collectors.joining("\t")));
//for each variant
stream().forEach(V->println(
//print the contig
V.getContig()+"\t"+
// print the position
V.getStart()+"\t"+
// print the REF alleles
V.getReference().getDisplayString()+"\t"+
// print the ALT alleles
V.getAlternateAlleles().stream().map(A->A.getDisplayString()).collect(Collectors.joining(","))+"\t"+
// loop over each type
all_types.stream().map(T->

String.valueOf(
// collect all genotypes
V.getGenotypes().stream().
//genotype must have type T
filter(G->G.getType().equals(T)).
count()) +
"\t"+
// collect all genotypes
V.getGenotypes().stream().
//genotype must have type T
filter(G->G.getType().equals(T)).
//convert to sample name
map(G->G.getSampleName()).
// join string
collect(Collectors.joining(","))
).
collect(Collectors.joining("\t")))
);
```

a one-liner:

```
$ java -jar dist/bioalcidaejdk.jar -e 'final List all_types = Arrays.asList(GenotypeType.values());println("CHROM\tPOS\tREF\tALT\t"+all_types.stream().map(T->"count."+T.name()+"\tsample."+T.name()).collect(Collectors.joining("\t")));stream().forEach(V->println( V.getContig()+"\t"+ V.getStart()+"\t"+ V.getReference().getDisplayString()+"\t"+ V.getAlternateAlleles().stream().map(A->A.getDisplayString()).collect(Collectors.joining(","))+"\t"+ all_types.stream().map(T-> String.valueOf( V.getGenotypes().stream(). filter(G->G.getType().equals(T)). count()) + "\t"+ V.getGenotypes().stream(). filter(G->G.getType().equals(T)). map(G->G.getSampleName()). collect(Collectors.joining(",")) ). collect(Collectors.joining("\t"))) );' src/test/resources/rotavirus_rf.vcf.gz

CHROM	POS	REF	ALT	count.NO_CALL	sample.NO_CALL	count.HOM_REF	sample.HOM_REF	count.HET	sample.HET	count.HOM_VAR	sample.HOM_VAR	count.UNAVAILABLEsample.UNAVAILABLE	count.MIXED	sample.MIXED
RF01	970	A	C	0	4	S1,S2,S3,S4	0	1	S5	0	0
RF02	251	A	T	0	3	S1,S4,S5	2	S2,S3	0	0	0
RF02	578	G	A	0	4	S1,S2,S3,S5	0	1	S4	0	0
RF02	877	T	A	0	4	S2,S3,S4,S5	1	S1	0	0	0
RF02	1726	T	G	0	3	S1,S4,S5	2	S2,S3	0	0	0
RF02	1962	TACA	TA	0	4	S2,S3,S4,S5	1	S1	0	0	0
RF03	1221	C	G	0	3	S1,S4,S5	0	2	S2,S3	0	0
RF03	1242	C	A	0	4	S1,S2,S3,S5	0	1	S4	0	0
RF03	1688	T	G	0	4	S1,S2,S3,S4	0	1	S5	0	0
```
Question: Hi all,

I am currently trying to perform a new assembly using paired and mate paired reads. I first mapped my paired-end on a reference genome. Now I would like to extract from the mapping only the covered portions and get contigs from my reads (not a consensus with the reference genome including N between the uncovered portions). I would like then to use the paired-end reads and the mate paired to scaffold my assembly.

I there a proper way to do that ? I think maybe mpileup could help, but until now I only managed to produced consensus.

Extra question, which software would you recommend to perform the scaffolding ? Does any one exist allowing different insert size ?

Thank you very much!

Clément

Answer: Here is an approach that might do what you want, I've tested it with some random data - In the first panel are mappings of some simulated reads with some SNPs, below are the mappings of the two generated contigs.

To generate the consensus contigs I used two modules from proovread, a correction software I have written. You will need to clone/install it from github and you will need a recent version of samtools for it to run.

# generate consensus (ignoring reference)
bam2cns --bam reads-to-ref.bam --ref ref.fa --prefix cns.fq

# trim no-coverage regions and convert to fasta
SeqFilter --in cns.fq --out cns.trimmed.fa \
--fasta --min-length 100 \
--phred-offset 33 --trim-win 1,1
Question: A quick question: I want to count the total aligned bases in a BAM file. I try to apply pysam (samtools) and pybedtools (bedtools) to do that. How can I make it?
Thanks.

**EDIT**:
I have solved this question with `pybedtools` and `pandas`:

genomecov_base = bd.BedTool.genome_coverage(bam, dz=True) ## Depth per base on genome

gencov_perbase = pan.read_table(genomecov_base.fn, names=['chrom', 'position', 'depth'])
gencov_allbase = (gencov_perbase.iloc[:, 2]).sum() ## Total bases mapped on genome

Answer: Hi, if it is available you could do a quick bedtools and awk calculation:

bedtools genomecov -ibam input.bam -bga > output_aln.bam
awk '{s+=$4}END{print s}' output_aln.bam

this will output the total of reads at each base from your reference genome in your bam file, and then awk will just some the number of reads at each base (giving you the total aligned bases).

I hope this helps!
Question: I am working with a program that generates artificial FASTQ files from a given reference genome and it allows for read length customization. I am trying to produce several FASTQ files with different read lengths. Hence, I want to know what read lengths are produced by real sequencers such as those from Illumina. I found [this][1] site which states the maximum read length of some sequencers.

However, when stating, for example, that the maximum read length is 2 x 150 bp, is that supposed to mean that each end of a contig will be 150 bp long? Also, if that is the case, does this mean that the sequencer could be set up to produce read lengths equal to any number of base pairs less than 150? i.e. could the sequencer produce read lengths of 2 x 1 bp, 2 x 2 bp, 2 x 3 bp,...., 2 x 148 bp, 2 x 149 bp? Lastly, which one of those sequencers is the most commonly used for genomes which are about as long as, say, that of E. Coli?

[1]: http://www.illumina.com/systems/sequencing.ilmn

Answer: *However, when stating, for example, that the maximum read length is 2 x 150 bp, is that supposed to mean that each end of a contig will be 150 bp long?*

No. Contigs are [something else][1] based on an assembled consensus sequence.

2 x 150 bp means that you have two 150 bp reads of sequence data from a single piece of DNA. The pair of reads are separated by an unspecified length based on the insert size (usually 200 to 1200 bp) which is size-selected during the sequencing prep. This picture might make things more clear.

[Mapping Reads][2] by [Suspencewl][3] - Own work. Licensed under [CC0][4] via [Wikimedia Commons][5].

The blown up read is 2 x 35. The insert size (bp between the sequencing adapters) is 400-500 bp. All of the reads in the picture could be used to assemble a single contig based on the consensus sequence.

For illumina sequencing, the read length is specified by the reagent kit, so you have limited flexibility there.

The MiSeq is capable of 15 Gb of sequence output. The *E. coli* genome is ~4.6 Mb in length so in a single run a MiSeq could easily sequence an *E. coli* genome at ultra-high depth.

[1]: http://genome.jgi-psf.org/help/scaffolds.html
[2]: http://commons.wikimedia.org/wiki/File:Mapping_Reads.png#mediaviewer/File:Mapping_Reads.png
[3]: http://commons.wikimedia.org/wiki/User:Suspencewl
[4]: http://creativecommons.org/publicdomain/zero/1.0/deed.en
[5]: http://commons.wikimedia.org/wiki/
Question: Hi,

I have mapped (paired-end) RNA-seq reads (so BAM file information) and I want to extract the nucleotide in the genome immediately prior to the (forward) mapping (essentially to test whether a RNAse digest has worked correctly).

As it is in BAM format, then extracting the position of the read (chromosome and location) is straightforward (just samtools view it, and process the line output for the relevant field), but I was wondering if there was a quicker way than brute-force (i.e. extracting, sorting, then reading from the genome FASTA file).

Answer: I think the easiest approach would be to just add an "N" to both ends of the reads before you map them. Then you can see the reference base in the MD tag.
Question: I am aiming to determine the species/subspecies identity of a (possibly bacteria) sample. I applied the following [sendsketch][1] command:

sendsketch.sh in=contigs.fasta,

where the input file (contigs.fasta) was previously generated from the raw paired-end read file using SPAdes. The results of the above command can be viewed here: [here][2]. The sendsketch website states that WKID is "the column that tells you the actual sequence similarity (disregarding size)." As a result, I believe I should focus on H. macacae since it had the highest WKID value (along with KID = 0.05%; ANI = 96.6%; Contam = 2.5%).

1) Are the output values here sound or a cause for alarm? Is the WKID (38.9%) or KID (0.05%) value too low? Or are other values questionable?

2) If the values are sound, I hope to check the mapping rate (using BWA) of the sample to the H. macacae genome (so that I can use it as a reference genome for further steps). However, I am a bit stuck on figuring out the appropriate sequence to download. The taxID number output from the sendsketch was 398626. It seems there are 46 sequences available on NCBI ([here][3]). How should I go about determining the optimal H. macacae sequence to map my samples to?

Side note: I have a similar problem with H. mastomyrinus, the second-highest WKID value from sendsketch. There are 17 sequences available on NCBI [here][4]. In this case, some are 16S and some are 23S. How do I know which sequence to map against in this scenario?

[1]: http://seqanswers.com/forums/showthread.php?t=74019
[2]: https://ibb.co/HpT4nJP
[3]: https://www.ncbi.nlm.nih.gov/nuccore/?term=txid398626%5BOrganism%5D
[4]: https://www.ncbi.nlm.nih.gov/nucleotide?term=txid287948[Organism]

Answer: There's only one H. macacae genome in my slightly dated local RefSeq Bacteria copy:

GCF_000507845.1	PRJNA224116	SAMN00690723	AZJI00000000.1	representative genome	1357400	398626	Helicobacter macacae MIT 99-5501	strain=MIT 99-5501	latest	Scaffold	Major	Full	2013/12/13	Heli_maca_MIT_99-5501_V1	Broad Institute	GCA_000507845.1	identical	ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/507/845/GCF_000507845.1_Heli_maca_MIT_99-5501_V1	assembly from type material
Question: Hello biostars,

I'm wondering if there is a way to find all genes sharing the same protein domain, e.g. the **WD repeat domain containing** occurs in many genes: BUB3, CDC20, DDB2, etc.

Is there a way I can retrieve the full list of genes having domain XX in common?

Thank you in advance!

Answer: Any time you want to produce a list of some description, try using BioMart first:http://www.ensembl.org/biomart/

Once you've selected which gene dataset you want to look at (e.g. Homo Sapiens (GRCh38.p5) from Ensembl Genes 84), you can go to **Filters** > **PROTEIN DOMAINS AND FAMILIES** > **Limit to genes with these family or domain IDs** and type in your domain.

You'll find BioMart comes in handy for a large number of tasks.
Question: I followed the post related to [PCA from vcf file][1]
followed by

library("SNPRelate")
vcf.fn<-"~/xxx/tmp.vcf"
snpgdsVCF2GDS(vcf.fn, "ccm.gds", method="biallelic.only")
genofile <- openfn.gds("ccm.gds")
ccm_pca<-snpgdsPCA(genofile)
plot(ccm_pca$eigenvect[,1],ccm_pca$eigenvect[,2] ,col=as.numeric(substr(ccm_pca$sample, 1,3) == 'CCM')+3, pch=2)

[1]: https://www.biostars.org/p/44735/#423942

I got a plot but not able to conclude from it. There is no legend or color code differentiation between the breeds. My understanding to R is very limited and couldn't able to find a fix. May I get some help/tips to fix this problem.
I will be grateful to you all.

Answer: I am confident that you know to what these numbers relate?

# set specific colours for each factor level

You can set colours manually to any vector like this:

x <- c('01','03','05','01','08','08','01','03','05')
unique(x)
[1] "01" "03" "05" "08"
colours <- c('royalblue', 'forestgreen', 'red4', 'pink')[factor(x, levels = c('01', '03', '05', '08'))]

Here, the order of the colours (`'royalblue', 'forestgreen', 'red4', 'pink'`) must match the order of the factor levels that you are setting (`'01', '03', '05', '08'`).

# apply a colour gradient / pre-built palette to your factor levels
You can also automatically assign colours via *RColorBrewer*:

require('RColorBrewer')
pick.col <- brewer.pal(10, 'Spectral')
colII <- colorRampPalette(pick.col)(length(unique(x)))[factor(x, levels = c('01', '03', '05', '08'))]

Here, I use the '*Spectral*' palette. For other palettes, take a look here (half-way): http://www.sthda.com/english/wiki/colors-in-r

# create your own palette

For a custom palette, do something like this:

myPalette <- c('white','pink','red2','red4')
colorRampPalette(myPalette)(length(unique(x)))[factor(x, levels = c('01', '03', '05', '08'))]

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

You can then plot a custom legend via the `legend()` function, which comes with base R.

Kevin