cannin/biostars_qa huggingface dataset:
Question: I have a little example fasta file with protein fragments to figure out how usearch works:

fragments

>seq1
TKEHALSKERAA

>seq2
KKEHALSKERAR

>seq3
AAHAASAERAAE

>seq4
AAHAASAERAAS

I used usearch with the following options:

usearch -cluster_fast ex.fasta -id 0.5 -uc cluster.uc

I expect at most 2 clusters (1,2) and (3,4) but the result contains only singletons.I decreased the identity and also the gap penalties but the result is the same.

Any idea?

Thank You

Answer: I think your sequences are too short. Try reducing the `-minhsp` value, but it simply may not work with these sequence lengths. If you duplicate each sequence and make the length 24 instead of 12, you will get the desired clusters.
Question: I have a very long list of coordinates that looks like that: `genome_build:chrN:coordinate`

I want to know if the given coordinate hits a gene or whatever else. What is the most efficient way to do so using python?

Example: having hg19:chr8:119331103 as input, it should return [SAMD12]

Answer: Hi,

You could use [Ensembl's rest API][1] for that using requests.

But if you have a very long list I would suggest that you try [pyGeno][2], it's typically the kind of applications if was made for:

```
x1 = 119331103
g = Genome(name = "GRCh37.75")
for gene in g.get(Gene, {"start >=": x1, "end <" : x1, "chromosome.number" : "8"}) :
print gene
```

[1]: http://http://rest.ensembl.org/documentation/info/overlap_id
[2]: http://pygeno.iric.ca/
Question: Hi there,

I have looked a bit around for this question but i haven't managed to fully understand the answers. I have done a WGCNA analysis on my data which ended up identifying several modules of co-expressed genes. Now, I would like to calculate significance of the correlation between the eigengenes and the trait data to further narrow down on what is more interesting. My traits data however has two covariates: time points, which is divided in 3, 7 and 35 days and virus strains, which are divided in 5 groups (control included). As far as i understand it shouldn't be allowed to replace integers for each strain group in order to calculate some kind of correlation, or is it? What other statistical analysis could i perform?

Answer: If these are the only traits in which you're interested, then you can either correlate the module values to these (with the traits encoded numerically as 1, 2, 3, etc), or, better, build a multinomial logistic regression model with the module values as x and time/strain as y ( `glm(time/strain ~ Module1)`; `glm(time/strain ~ Module2)`; *et cetera*)

You may also have to consider dividing up the analysis into multiple analyses, and contrast/compare the results manually. For example, running WGCNA separately for the different time-points may be an idea, and then building the regression model predicting for virus strain each time.
Question: I have few motifs whose location on a given DNA sequence has been identified using MATCH (from TRANSFAC). MATCH provides the location of these matrices in a numerical way as motif 1 236 (+) (for an example - motif 1 is present at location 236). I want to get a graphical output by feeding the length of DNA sequence and the the position of matrices so that its easy for presentation and comparing two different sequences with almost similar motifs

Thanks in advance for any suggestions

Answer: I always liked Melina II for graphical display of the features on the sequence, and it gives a sequence logo as well. http://melina2.hgc.jp/public/index.html But if your motifs aren't identified with the selected tools they offer for "methods" I don't know if it would work for you.
Question: I have downloaded a database from GEO using getGEO and then converted it into ExpressionSets using GDS2eSet. Now I am trying to normalize it using rma. It gives me error :unable to find an inherited method for function ‘probeNames’ for signature ‘"ExpressionSet". Which function I showed use to normalize the expression set ?

a1<-getGEO('GDS4758', destdir=".")
eset<- GDS2eSet(a1, do.log2=TRUE)
eset1<-rma(eset)

Answer: Hi,

the error message resembles one in this [thread][1]. It looks like you should use *different* **rma** function, namely from **oligo** package. You could try `oligo::rma` instead of `rma` (after possible installation of `oligo` library)

[1]: https://support.bioconductor.org/p/71097/
Question: I have a VCF file and I want to remove sites that are not polymorphic sites i.e. sites where all individuals are REF/REF or ALT/ALT. Is there a quick way to do this using bcftools or vcftools?

Thanks.

Answer: Now I hopefully understand your goal. This will remove side where all samples are either hom ref or hom alt:

$ bcftools view -e 'COUNT(GT="AA")=N_SAMPLES || COUNT(GT="RR")=N_SAMPLES' input.vcf

Have a look at the [manual][1] to explore more ways to filter `vcf` files using `bcftools`.

I've deleted my old answer, because it hasn't answer your question.

[1]: http://www.htslib.org/doc/bcftools.html#expressions
Question: I have sequencing data of five bacteria, which were generated using Illumina MiSeq. Four of them were sequenced using a paired-end 2x300 protocol and one was sequenced using the nextera mate-pair protocol.

My question is: What are the softwares that you recommend me to assemble these genomes (the largest has almost 8 Mbp)?

I have access to a CLC Workbench. It seems quite ease to use, but I dont know if it is the best one. Most of papers that I found that evaluate the performance of assemblers from two-three years ago.

I also have to mention that I have two i7 with 8 threads PCs available for this objective (one with 32 and another with 8 Gb RAM).

Thanks in advance.

Answer: You should read in the literature to know which one is the best one for you specific data. Here you have a nice paper comparing some assembly tools, and it is a recent paper (2014): http://genomebiology.com/2014/15/3/R42

In my opinion SOAPdenovo2 and SGA are a good choice. Bambus is quite difficult to install and to understand. SPACE also is nice, but if you want to use the last version you have to pay so...

Hope it helps.
Question: Hi

I have to run this on 136 .vcf files but manually one by one takes an ages

Can you please help me in any script to do that at the same time in terminal?

[fi1d18@cyan01 snp]$ bcftools query -f '%CHROM\t%POS\t%REF\t%ALT[\t%ID]\n' file.vcf > file.txt
[fi1d18@cyan01 snp]$

Thanks

Answer: Thank you @genomax

Similar to what you suggested, in terminal one should type

for f in *.vcf; bcftools query -f '%CHROM\t%POS\t%REF\t%ALT[\t%ID]\n' "$f" > ${f%.*}.txt done
Question: I have a gff file that contains all the information from Chromosome 1 - Chromosome 14. But I need gff information on individual chromosome basis. For example, I am performing some experiment on Chromosome 11 and trying to visualize my result on IGV. When I load the gff file in IGV it is showing gene information of all the chromosome. How can I get the gff of only Chromosome 11?

Answer: Assuming your chromosomes are named `chrNN` the following should extract `chr11`

grep chr11 your_file.gff > chr11.gff

If chr11 is in the first column and you want only those lines then do

grep ^chr11 your_file.gff > chr11.gff
Question: Hi All,

Do you have R code for the bisulfite sequencing PCR result plotting script.

Thanks

[][1]

[1]: http://www.springerimages.com/img/Images/BMC/VOL=2011.3/ISU=1/ART=4/MediaObjects/MEDIUM_13148_2011_4_Fig4_HTML.jpg

Answer: There might be a package that will produce those exactly, but if not you can do it in ggplot2:

```r
d <- data.frame(CpG=rep(1:10,4), read=rep(1:4, each=10), value=sample(c(1,16), 40, replace=T))
ggplot() + scale_shape_identity() + geom_point(data=d, aes(x=CpG, y=read, shape=value), size=5)
```

Note that you'll need to tweak the size, remove the background, etc., but that should serve as a start.
Question: Hello everyone,

I was thinking, if for example I can use samples generated from **Illumina HumanHT-12 V4.0 expression beadchip** as a healthy control for samples generated in **Illumina HumanHT-12 V3.0 expression beadchip** or vice versa.

Is this possible or not? if the answer is no, what are the requirements to use data generated in other microarray experiments as a control??

Best regards,

Leite

Answer: The short answer is no. For the longer answer, you'll need to go with the age old rule of thumb, keep all experimental conditions the same, except the variable that's relevant to your hypothesis. In the realm of sequencing, and high throughput experimentation, technical variability plays a huge role.

Three things to consider when designing a sequencing / high throughput experiment: batches, platform, and sample type placement.

To absorb a given *nuisance* effect, there needs to be balance on your primary term of interest. Take the following design as an example :

> pheno
Sample ID SampleType
SAM1 A
SAM2 A
SAM3 A
SAM4 B
SAM5 B
SAM6 B

These are all from the same sequencer, no nuisance variables that need to be accounted for, so the design matrix is simply `design.matrix(~0 + SampleType, data = pheno)`. Now consider a known batch effect such as this:

> pheno.batch
Sample ID SampleType Batch
SAM1 A B1
SAM2 A B1
SAM3 A B2
SAM4 A B2
SAM5 B B1
SAM6 B B1
SAM7 B B2
SAM8 B B2

We'd consider this experiment to be well balanced, as you can see there are samples with `SampleType` `A` and `B` in batch 1 and 2. This balance means that variation can be estimated across `SampleType` and `Batch`, with the following design matrix: `design.matrix(~0 + SampleType + Batch, data = pheno.batch)`.

Next, here's an example of an unbalanced design, where the design matrix will not be `full rank`, which means that there's an unbalanced term(s) in the design matrix. In this case, we can see that our `SampleType` column, and `Batch` columns are the same when looking at factor levels.

> pheno.batch2
Sample ID SampleType Batch
SAM1 A B1
SAM2 A B1
SAM3 A B1
SAM4 A B1
SAM5 B B2
SAM6 B B2
SAM7 B B2
SAM8 B B2

When running `design.matrix(~0 + SampleType + Batch, data = pheno.batch2)`, there will be a `full rank` error. Now that we understand the difference between a balanced and unbalanced design, lets take your OP as an example:

> pheno.batch.op
Sample ID SampleType Platform
SAM1 Healthy_Control HT12v4
SAM2 Healthy_Control HT12v4
SAM3 Healthy_Control HT12v4
SAM4 Healthy_Control HT12v4
SAM5 Disease HT12v3
SAM6 Disease HT12v3
SAM7 Disease HT12v3
SAM8 Disease HT12v3

In this case, `Platform` is essentially, the same as `SampleType`, so we can't untie that variation. What that essentially means at an abstract level, is that we can't statistically tell the difference between variation that is coming from the difference of platform, or biological variance coming from the sample types. These experiments are highly sensitive to change, so that's not even accounting for library preps, sample preparations, source biological material, kits used, temperature in the room, person who prepped the sample, etc.
Question: Hi all,

Currently, I am working on my very first RNAseq study and have met a dilemma where inputs from more experienced bioinformaticians would be amazing.

For a differential gene expression study in a non-model organism, a de novo reference transcriptome was assembled from 300 M reads in Trinity. For 3 experimental conditions (1 negative control, 1 positive control and the treatment of interest) triplicate samples were sequenced with a depth of 25 M reads.

The reference transcriptome was annotated with Trinotate.

For differential gene expression determination the Kallisto/Sleuth pipeline was being used - and here comes my dilemma of best practices:

A number of the Trinity transcripts could not be annotated by Trinotate (NA) and is being dropped in the Sleuth analysis when using the *"so <- sleuth_prep(s2c, ~treat, target_mapping = annotation, aggregation_column = 'gene')"* expression.

I played around with the annotation file and replaced the NA's in the gene column with the corresponding Trinity transcript IDs, which included some of the transcripts as significantly differentially expressed.

**What is the right thing to do?**

- Would you let Sleuth drop the non-annotated transcripts, even though some of them are significantly differentially expressed?

- Or, would you include these transcripts in the "gene" column with their corresponding Trinity transcript IDs, even though they cannot be analyzed on the gene-level (the transcript isoforms cannot be collapsed in the analysis like the annotated ones)?

Thank you!

Cheers, Birgitte

Answer: In addition to what [andrew.j.skelton73][1] said, read carefully Trinity and Trinotate wikis, they have suggestions on how to filter RNAseq assemblies for further analyses. For differential expression, you should filter by level of expression, removing lowly expressed transcripts regardless of annotation status, as they will have low power anyway - specially with only three replicates per treatment.

[1]: https://www.biostars.org/u/9912/
Question: Hi guys,

I have R programming question: I have more than 1000 genes names `Names`. I want to match the samples containing those genes in `Names` and paste the corresponding gene with the ".AD" extension next to it in the same order to get the `Result` as shown below. Thank you.

```r
Names <- c("cebi", "pithe", "MAPK", "sapiens", "JUNK", "calli", "STR")
samples <- c("MAPK", "JUNK", "STR")
```

`Result`:

"cebi", "pithe", "MAPK", "MAPK.AD", "sapiens", "JUNK", "JUNK.AD", "calli", "STR", "STR.AD"

Answer: [EDIT: you clearly said that order matters; I didn't check back after posting into R. I'll leave this up in case someone can find some utility in it since another has provided an answer that works.]

Wasn't clear if order matters, but here's what you want:

```r
Names <- c("cebi", "pithe", "MAPK", "sapiens", "JUNK", "calli", "STR")
samples <- c("MAPK", "JUNK", "STR")
names <- Names
samples %in% names
[1] TRUE TRUE TRUE

match(samples, names)
[1] 3 5 7

names[match(samples, names)]
[1] "MAPK" "JUNK" "STR"
paste(names[match(samples, names)], ".AD", sep="")
[1] "MAPK.AD" "JUNK.AD" "STR.AD"
c(paste(names[match(samples, names)], ".AD", sep=""), names)
[1] "MAPK.AD" "JUNK.AD" "STR.AD" "cebi" "pithe" "MAPK" "sapiens"
[8] "JUNK" "calli" "STR"
sort(c(paste(names[match(samples, names)], ".AD", sep=""), names))
[1] "calli" "cebi" "JUNK" "JUNK.AD" "MAPK" "MAPK.AD" "pithe"
[8] "sapiens" "STR" "STR.AD"​
```
Question: I have used `featureCounts` to get the read counts one a gene level with the following command:

featureCounts -a gencode.gtf -F GTF -p -s 2 -T 8 -o sample.txt sample.sorted.bam

I'm interested in extracting the counts on a transcript level. My data is paired end and reverse stranded. I'm a bit confused with command I have to use to get the counts on transcript level [featureCounts manual][1]

Can anyone please tell me how to get the counts on transcript level? thanq

[1]: http://bioinf.wehi.edu.au/featureCounts/

Answer: You'll probably want to use `RSEM` or `stringTie`
Question: Hello,

I would be grateful for advice on good ticketing systems for bioinfo support. Do standard IT ticket systems like Fork8 etc are appropriate to use in genomics or. Earlier I used gitlab issues but the clients also needed to have a gitlab account. I can think of a few requirements but please feel free to add to them.

- Clients can no see each others tickets (this excludes deploying biostars instance)

- Good formatting support for code and comments

- Possible git/gitlab integration

Also if the ticketing system has helped you in your evaluation/tenure/promotion please let us know.

Answer: We use [Request Tracker](https://bestpractical.com/request-tracker) in Ensembl. All tickets are private and can be conducted entirely by email, so no login is required by requestors. You can also comment on a ticket without a requestor seeing it, which can be useful for keeping track of your own internal investigations.

In RT, you "take" a ticket, then you are the requestor's only point of contact. You can change the settings so that everybody who follows the queue sees all emails, or so that only the ticket owner sees all subsequent emails. You can "steal" tickets from someone if it becomes apparent that someone else is better equipped to deal with it, or if the owner goes on holiday or something.

It's not great for formatting. Generally the interface is a bit clunky. You can format your messages to include code and stuff, but it's certainly not as easy to use as something like BioStars. Since the requestor is using email, you suffer from all kinds of dodgy formatting that their email client will do to you. Mostly I send stuff like scripts as attachments.

I'm not sure if you can integrate with git/gitlab. It has an API to allow you to write custom scripts in Perl to do quite a lot of stuff, so it's probably possible.

It's fairly easy to search and generate graphs and stats of your search results. I personally haven't used RT to help justify a promotion or similar, but you could certainly show that an individual took some percentage of tickets and resolved them all within a certain space of time (although this tells you nothing about the difficulty of those tickets, it could be like comparing the number of operations and survival rate of a brain surgeon to one who mainly does C-sections).
Question: Dear all,
I want to ask a very basic question. I am making a bubble plot using ggplot2 having table structure:

GO term Number Class Type
1. Metabolism 5 start duf BP
2. Photosynthesis 10 hzs BP
3. Nucleus 15 hs CC
4. Kinase 16 hs MF

I want to make a bubble plot having Number on x axis, GO term on y axis, bubble color should be based on Class and background color will be based on Type.
My R code is:

> ggplot(bubble_plot, aes(x=Number, y=GO term, size = Number, col =
> Class)) + geom_point(alpha=0.7)

With this I am getting the desired plot except the background color, When I am using "fill = Type", I am not getting the background color based on the BP, CC or MF based on the " Type" column in the table.

The desired plot should look like:

Please help.

Thanks in advance

Answer: So with your data that seems to look something like this:

structure(list(GO_term = structure(c(2L, 4L, 3L, 1L), .Label = c("Kinase",
"Metabolism", "Nucleus", "Photosynthesis"), class = "factor"),
Number = c(5L, 10L, 15L, 16L), Class = structure(c(3L, 2L,
1L, 1L), .Label = c("hs", "hzs", "start_duf"), class = "factor"),
Type = structure(c(1L, 1L, 2L, 3L), .Label = c("BP", "CC",
"MF"), class = "factor")), class = "data.frame", row.names = c(NA,
-4L))

You could work with `geom_tile` and set its `width = Inf` to get something akin to that plot you're trying to emulate. The use of `forcats::reorder()` within `aes()` is important in order to group the Y-axis values together on the basis of your `Type` column.

library(ggplot2)
library(forcats)

ggplot(mydat, aes(y = reorder(GO_term, as.numeric(Type)), x = Number, size = Number)) + geom_point(aes(color = Class), alpha = 1.0) +
geom_tile(aes(width = Inf, fill = Type), alpha = 0.4) +
scale_fill_manual(values = c("green", "red", "blue"))

Which yields:
![col_by_bg2.png][1]

The problem you'd probably run into is having to pass the appropriate number of colors to `scale_fill_manual` for your actual dataset.

[1]: https://raw.githubusercontent.com/vragh/share/master/col_by_bg2.png
Question: Hi!!

I have to analyze paired-end RNA-seq read that are in an unusual format: both pair-end reads are joined in one FASTQ. I need to split the file in two separated FASTQ paried-end files.

There are a galaxy tool named FASTQ splitter that can do this:

---

FASTQ splitter

What it does

Splits a single fastq dataset representing paired-end run into two datasets (one for each end). This tool works only for datasets where both ends have the same length.

Sequence identifiers will have /1 or /2 appended for the split left-hand and right-hand reads, respectively.

Input format

A multiple-fastq file, for example:

@HWI-EAS91_1_30788AAXX:7:21:1542:1758
GTCAATTGTACTGGTCAATACTAAAAGAATAGGATCGCTCCTAGCATCTGGAGTCTCTATCACCTGAGCCCA
+HWI-EAS91_1_30788AAXX:7:21:1542:1758
hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh`hfhhVZSWehR

Outputs

Left-hand Read:

@HWI-EAS91_1_30788AAXX:7:21:1542:1758/1
GTCAATTGTACTGGTCAATACTAAAAGAATAGGATC
+HWI-EAS91_1_30788AAXX:7:21:1542:1758/1
hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh

Right-hand Read:

@HWI-EAS91_1_30788AAXX:7:21:1542:1758/2
GCTCCTAGCATCTGGAGTCTCTATCACCTGAGCCCA
+HWI-EAS91_1_30788AAXX:7:21:1542:1758/2
hhhhhhhhhhhhhhhhhhhhhhhh`hfhhVZSWehR

---

Do you know any other standard alone script that can do this job?

Answer: copy this code to a file, e.g. `splitPairedEndReads.pl`

use strict;
use warnings;
my $file = $ARGV[0];
open(FILE, "<$file") || die "cannot open $file\n";
open(OUT1, ">$file\_1") || die "cannot open $file\_1\n";
open(OUT2, ">$file\_2") || die "cannot open $file\_2\n";
while(){
chomp;
print OUT1 "$_\/1\n";
print OUT2 "$_\/2\n";
my $newline = ; chomp($newline);
print OUT1 substr($newline, 0, length($newline)/2)."\n";
print OUT2 substr($newline, length($newline)/2, length($newline)/2)."\n";
$newline = ; chomp($newline);
print OUT1 "$newline\/1\n";
print OUT2 "$newline\/2\n";
$newline = ; chomp($newline);
print OUT1 substr($newline, 0, length($newline)/2)."\n";
print OUT2 substr($newline, length($newline)/2, length($newline)/2)."\n";
}
close(FILE);

Now call it with `./splitPairedEndReads.pl input.fastq`.

I hope it works.... haven't checked it.
Question: Hello all,

I have a GWAS cohort and a WGS cohort that is a subset of the GWAS cohort.

I have identified all the lead SNPs from the GWAS data (larger n) and have defined linkage disequilibrium around each lead SNP in the study.

Further, I have phased the WGS genomes. As such, I can determine the variants that are in phase, or not in phase, with the lead SNPs in the study.

Right now, it is my intention to do this by parsing the VCF files, then reformatting into SKAT-O format. However, my question is, *are there existing workflows that port phased genomes directly into any rare variant association tool, or, is one or more of these steps DIY?*

For clarity, what I mean by rare variant association tools is roughly summarized by [the Lee / Lin review paper in AJHG][1]

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

Answer: I wrote a script to do exactly this:

https://gist.github.com/brentp/4ec866b5ab3dd4e7bc46

It parses the VCF in python and sends the haplotype blocks to R using [pyper][1]. It can be a start for the parsing/creating haplotypes or you can use it as-is to get p-values from SKAT.

There are a lot of methods for testing haplotypes, but I think [this poster][2] does a nice job comparing.

[1]: https://pypi.python.org/pypi/PypeR/1.1.0
[2]: http://srcos2014.blogs.rice.edu/files/2014/07/AnandaDatta.pdf
Question: Hi Biostars,

I am looking for a tool to introduce variants in a reference genome.

**The ideal tool I am looking for:**

*Takes as input:*

- genome fasta

- a vcf file

*What it does:*

- Insert some of the variants (randomly) present in the vcf into the reference genome

*What is the output*

- The ideal output would be a mutated new reference fasta, and a vcf file (The vcf file is not that important).

The reason I am looking for this is because I am working on simulation tool, but I would be happy to rely on a third party software to introduce the mutations.

cheers,

Answer: Not exactly what you are asking but if you are not wedded to specific vcf files then you could go the other way by using `mutate.sh` from BBMap suite. Lots of processing parameters you can control.

$ mutate.sh

Written by Brian Bushnell
Last modified December 13, 2018

Description: Creates a mutant version of a genome.

Usage: mutate.sh in= out= id=
I/O parameters:
in= Input genome.
out= Output mutant genome.
vcf= VCF file showing variants added.