cannin/biostars_qa huggingface dataset:
Question: Hi guys, do you know based on which database does prokka assign EC numbers to found proteins? KEGG or NCBI database or an internal db?
thanks

Answer: Prokka assigns all annotations by copying fasta headers from its internal database, or from HMM descriptions.

The meaning of EC numbers is explained at [**ExPASy**][1], specifically in [**this file**][2]. Prokka has all the files used for annotations in `db/genus` and `db/kingdom` directories.

[1]: https://enzyme.expasy.org/
[2]: ftp://ftp.expasy.org/databases/enzyme/enzclass.txt
Question: Hi,

I believe a similar question has been asked before ([Question: Is Star Managing Strand-Specific Data ?
][1]), but in the answer to that, it was mentioned that we don't need to do anything extra in STAR alignment if we have strand-aware RNA-seq data. But at another place ([Are any additional parameters for STAR to hand the stranded RNA-Seq data?][2]), the author of STAR mentions that presently STAR doesn't use strand information for mapping, and that if we use a certain option (--quantMode), it provides output for 3 kinds of strandedness).

So which is it? Does STAR perform strand-aware mapping, which we can use downstream to get antisense transcripts, or if not, then what is the best way to achieve that? I am working with Tru-Seq Stranded Total RNA library, which apparently maps to the opposite strand of the genome (opposite to what you would normally expect).

I am primarily interested in finding novel non-coding transcripts in my data

[1]: https://www.biostars.org/p/62974/
[2]: https://groups.google.com/d/msg/rna-star/oGvChGG2gto/UQrP3c8yDQAJ

Answer: >strand-aware RNA-seq data

There is no such thing. There are strand-specific library preparations, though, which generate strand-specific (relative to the mRNA) reads, which is probably what you mean.

> Does STAR perform strand-aware mapping, which we can use downstream to get antisense transcripts

You can't get transcripts, sense or antisense - you can get mappings and counts from STAR. What are the downstream analyses you intend to perform?

STAR maps reads to an indexed (most commonly annotated) reference genome, and it doesn't care about transcript strandedness. From these mappings, it can output many different things, which can be used for several downstream analyses. The three most commonly used outputs from STAR are:

1. read mappings to the genome, in BAM format. By default, strandedness has no importance: STAR maps the reads to the best location, regardless of strand, gene location or whatever. However, the parameter `--outSAMstrandField intronMotif` - **which should be used only if you want to analyse downstream with Cufflinks, and have unstranded RNAseq dara** - alters the default behaviour and adds XS strand attribute. Read the section **Compatibility with Cufflinks/Cuffdiff** from STAR's manual. If you have stranded data, **you don't have to change anything in STAR**.

2. read mappings to the transcriptome, in BAM format. Those are to be used with RSEM, mostly.

3. counts of reads mapping to genes, similar to HTSeq or featureCounts. Here STAR outputs three columns of counts, one for unstranded, one for stranded, and one for reverse stranded.

Most of these points were already clearly made by Alex Dobin in the thread [Are any additional parameters for STAR to hand the stranded RNA-Seq data?][1], and the STAR manual.

[1]: https://groups.google.com/forum/#!msg/rna-star/oGvChGG2gto/UQrP3c8yDQAJ
Question: Hi all,

I have a genome-wide list of germline SNPs and short indels for Arabidopsis thaliana, which I generated with Varscan. Regardless of the tool used to generate them, I would like to annotate them, i.e. knowing which ones can cause an aminoacid change or an early stop, using the default Arabidopsis thaliana Columbia 0 cultivar, for which I have both the sequence (FASTA from TAIR10) and the updated annotation (GFF from Araport).

Chrom	Position	Ref	Var
Chr1	626503	G	T
Chr1	926694	C	T
Chr1	5280350	C	A
Chr1	5699993	C	A
Chr1	7004559	G	A
Chr1	8325810	C	T
Chr1	9371723	T	G

What I want to do is similar to what Annovar does, but unfortunately Annovar does not support Arabidopsis. I was thinking of an already existing R pipeline that takes in a genome, an annotation, a SNP/indel list and *boom*, annotation. But I couldn't find any, except maybe [snpEffect][1]. Any tips? Thanks in advance!

[1]: https://github.com/douglasgscofield/snpEffect

Answer: Hello,

have a look at [snpEff][1]. Your list of SNP/indel have to be in vcf or bed format.

fin swimmer

[1]: http://snpeff.sourceforge.net/SnpEff_manual.html
Question: Dear All:
Does anyone know an efficient tool or script that removes certain reads with from Fastq file when it matches specific sequence?
this means : if the reads contain my sequence so the entire read should be completely removed from the FASTQ file.

Many many thanks

Answer: Use [seqkit grep][1]:

$ zcat input.fq.gz | seqkit grep -v -s -i -p aggcg

This will remove all reads containing "aggcg".

[1]: https://bioinf.shenwei.me/seqkit/usage/#grep
Question: Hi everybody,

I am currently working on several projects using Illumina 450k DNA Methylation Microarrays. In order to detect Differentially Methylated Probes (DMP), I usually employ the Empirical Bayes-based method in the limma package. Based on the paper by Pan Du et al, I stick with M-values to fit the statistical model, and usually employ beta-values just in graphs and reports that are going to be read by fellow biologists.

In order to retain only results with a certain biological relevance, we usually apply a threshold on effect size, keeping only the probes that are significant and with a big effect size. Nothing new or fancy to this point. However, it is not uncommon to argue with fellows at the lab about the benefits or drawbacks of M-values against beta-values when trying to set a coherent effect size threshold.

Mathematical properties of M-values let us set a fixed threshold on effect size, while this is not that easy for betas. However, setting a threshold on M-values differences (say 1.4, as stated in the previous paper) usually results in a set of probes where a lot of them seem to present very small differences in beta-values (say, for example, 0.01), specially those near the minimum and maximum. This is very counterintuitive for a biologist, who is going to argue against it based on the fact that such a small difference means nothing from a biological point of view.

My mental idea of what's happening there is related to the technical bias of the 450k. This is, I try to convince people that a small difference in that region is as credible as a bigger difference in the middle region (around 0.5 in beta), due to the array design, but I do not think if I am right, or if I can even see correctly what is going on.

What do you usually do in your pipelines? Use M-values for the fit, and beta-values differences as thresholds for effect sizes? A first threshold on M-values and a second filtering step based on betas? Everything with betas? Nothing at all?

Any hint will be much appreciated.

Answer: Using M-values for statistics is by far the best way to go.

Regarding the utility of beta values, it's good to filter by those as well. Why? Because a 1% change in methylation is unlikely to be biologically meaningful. I've seen a number of papers publish such changes, but it's always interesting to note that they never bother to show functional relevance (probably because there is none). If you can show the biological relevance of such small changes then go ahead and follow up on them. Personally, I want to prioritize the results on the likelihood that it's causing some phenotype and using methylation changes for that (and only that) makes sense.

Edit: There's corollary to other types of data. For example, we can often detect small changes in highly-expressed genes when we do RNAseq. Some of these are meaningful, but the really small changes aren't. I wouldn't toss these results, but I also wouldn't recommend anyone pursue them first when doing a follow-up.
Question: I would like to get the r-squared for two SNPs in my population. One of the SNPs is in a mySQL db, and the other is currently also in a mySQL db, but I have the plink files for that one.

Can anyone advise me on how I will go about getting these into a format that Haploview will understand?

EDIT:

Just to clarify, I've two SNPs from different data sources on the patients. While they are indicated to be perfect proxies for each other (1000Genomes), they present different p-values when the Jonckheere-Tempstra test is performed, SNP-A (p=0.2192), SNP-B (p=0.0152). I've been asked to get the r-square between them.

Haploview accepts a plink format, I don't have the plink files for SNP-A, that is contained entirely within a mySQLdb, converted over from the SAS format of a previous student. I do have the plink files for SNP-B, but it is split into two separate plink datasets. Both datasets can be merged using plink, but there is still the issue of SNP-A not being available in a plink format.

Answer: If you have genotype for both SNPs for the same samples, simply create [PED MAP][1] files, (using ~~Excel~~, text editor, command line...) then use plink for [merging][2] and [calculating LD][3].

[1]: http://pngu.mgh.harvard.edu/~purcell/plink/data.shtml#ped
[2]: http://pngu.mgh.harvard.edu/~purcell/plink/dataman.shtml#merge
[3]: http://pngu.mgh.harvard.edu/~purcell/plink/ld.shtml
Question: Hi I'm studying about sequencing data analysis. I performed a series of analysis and i finally get multi-sample VCF file. This VCF file is 99 test group SNPs and 1 control group SNPs mixed. I want only test group SNP, so i tried to many select variants method. but its still remain control SNPs. Is there any methods to know the separate control SNPs and test group SNPs? I will be happy if anybody suggest me idea regarding this. Thank you.

Answer: I understand you want to separate a large multisample VCF file into 1 column control and 99 columns with Test group SNVs ?

Can't vcflib do this for you ?

https://github.com/vcflib/vcflib#vcflib

vcfremovesamples or vcfkeepsamples.
Question: As per the title: why?

Does the preliminary PCA step merely maximise the chances of identifying the clusters that tSNE / UMAP later identify?

When doing this, in tutorials, it seems that people blindly choose a number of PCs as input to tSNE / UMAP, without checking the amount of variation explained by the PCs. For example, 20 PCs may only account for 30% overall variation.

Assuming one follows a standard procedure for, e.g., scRNA-seq, is this one step too many. A pipeline could be:

1. normalisation of raw counts
2. transformation of normalised counts
3. PCA-transformation (may include an additional pre-scaling step)
4. tSNE/UMAP

Kevin

Answer: The primary reason for this given in the tSNE literature is for the sake of computational efficiency. See, for example, [vanderMaaten 2008](https://lvdmaaten.github.io/publications/papers/JMLR_2008.pdf). In essence, tSNE requires pairwise comparison of datapoints, so it can be incredibly computationally taxing on scRNA-seq datasets unless the dimensionality undergoes an initial reduction.

I can't speak to UMAP, I'm not familiar enough with its inner-workings, but I presume the initial PCA is done for similar reasons.
Question: Hey everyone!

I have a dataframe in R, the columns being Gene names and the rows being Isolate names. The value in each column can be 1 or 0, 1 for "Gene is present in Isolate" and 0 for "Gene not present in Isolate". So all in all, the dataframe shows in one-zero fashion which Isolate has which genes.

Now I want to do a network analysis to see which genes are most likely to co-occur, assess the strength of their connection and so on.

In R, I have tried the following:

```
>library(igraph)
>library(network)
>library(sna)
>library(ndtv)

>Genematrix <- data.matrix(df)
>g <- network(Genematrix, directed=FALSE)
> summary(g)
>plot(g)
```

What I get from this is a network object with 236 vertices. But what I actually want is the Genes as vertices ( 21 columns), so I can see the clusters and connections between them. In many tutorials I have seen that I need and edge list and a node list. the edge list is I think what I get from >g <- network(Genematrix, directed=FALSE), but I don't know how to get the node list.

Can anyone explain how to solve my problem and what I have to do to get the network I want?

Answer: Cols are the genenames and you want a gene-gene adjacency matrix to feed into igraph/network etc. You could do this with an incidence matrix (rows=vertices, cols = edges) as well, and your current dataframe looks almost like an incidence matrix (but isn't, since some columns may have more than 2 entries, so don't represent edges; and indeed, the vertices you want in your graph are present in the columns)

The following should convert your dataframe into an adjacency matrix (the edge between vertices u and v being weighted by the number of isolates where they were both 1)

m <- as.matrix(df) # note that df is a base function in R, so isn't a very good variable name
adj.m <- t(m) %*% m

Note that `adj.m` has diagonal entries that are of no value in your analysis so I think you can do something like the following to get rid of them

diag(adj.m) <- 0

Russ
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: You can use bioMart, these two links may help you out

http://www.ensembl.info/blog/2015/06/01/biomart-or-how-to-access-the-ensembl-data-from-r/

The snp content is on page 12 of the bellow link

http://www.bioconductor.org/packages//2.7/bioc/vignettes/biomaRt/inst/doc/biomaRt.pdf

Or you can also go to http://www.ensembl.org/biomart/martview/b610c4828014c048f6bb2e1b651ce2d2 the on-line tool.
Question: Hi, this may sounds a bit trivial question, i indexed my ref genome, along with the indexed genome, there is also a file called .fasta.fai. My question is what is in the fasta.fai? I opened this file, it does not have any header. Did I do anything wrong ?

I then used the ref genome to call SNPs. I assume some sort of position information is used to tell me the position of snps in the ref genome. Which column is position? Thanks.

Here is a subset of the fasta.fai file

gi|394055774|gb|AGTA02000001.1| 5516 93 70 71
gi|394055773|gb|AGTA02000002.1| 2292 5781 70 71
gi|394055772|gb|AGTA02000003.1| 4668 8199 70 71
gi|394055771|gb|AGTA02000004.1| 1190 13027 70 71

Answer: The fasta.fai is the fasta index, and the one you posted looks legit.

For each row:

Column 1: The contig name. In your FASTA file, this is preceeded by '>'

Column 2: The number of bases in the contig

Column 3: The byte index of the file where the contig sequence begins. (Notice how it constantly increases by roughly the amount in column 2?)

Column 4: bases per line in the FASTA file

Column 5: bytes per line in the FASTA file
Question: I follow the instruction from https://www.bioconductor.org/install/. But I got error as following

> source("http://bioconductor.org/biocLite.R")
Warning in install.packages("BiocInstaller", repos = a["BioCsoft", "URL"]) :
'lib = "/share/apps/R/3.1.1/intel/2013.0.028/lib64/R/library"' is not writable
Error in install.packages("BiocInstaller", repos = a["BioCsoft", "URL"]) :
unable to install packages
Calls: source ... eval.parent -> eval -> eval -> eval -> eval -> install.packages
Execution halted
Can someone help me here?

Answer: `'lib = "/share/apps/R/3.1.1/intel/2013.0.028/lib64/R/library"' is not writable`

Obviously, you don't have the permission to install package to the R library folder.

1. add the following lines to `~/.Renviron`: `R_LIBS_USER='~/R/library'`

2. try again:

```
source("https://bioconductor.org/biocLite.R")
biocLite("snpStats")
```

You should be able to install any package you want.
Question: I have SNP IDs in this format: `rs3094315, rs3115860`.

Now, I want to map these snp ids to gene symbols in which this snp is present.

Can some one help me in this?

Answer: Probably this may answer your question:

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

Else, just have a look at this blog and also comments where they have talked about various methods.

http://www.gettinggeneticsdone.com/2011/06/mapping-snps-to-genes-for-gwas.html

If you SNP ids aren't huge, you can try some web based tools like: http://pfs.nus.edu.sg

Or use the [perl modules][1] or downloads from UCSC browser (explained in the blog and also in the biostar link I have provided)

[1]: https://github.com/inti/FORGE/wiki/Tutorials-1%3A-Make-snp-to-gene-mapping
Question: Dear friends,

I am very new to next gen sequencing. Using mplieup and bcftools in samtools I have generated a vcf file. I want to write a perl script to extract only the heterozygous snp. How can I understand which are heterozygous snps and how they can be extracted through a perl or any programming script.

I will be very grateful if anybody can shed some lights towards the doubt please.

Regards,

Shubhankar

Answer: GT tag with 0/1 is Heterozygous. You dont need to even write a script for that.

vcftools --gzvcf file.vcf.gz --extract-FORMAT-info GT | grep "0/1"

will work
Question: Hello,
I am having a hard time opening a very large bedgraph. I have been suggested to split my bam file by chromosome with chrom-bed.py but it didn't work.
Is there any other alternative?
Thanks,
GP.

Answer: I wrote a java tool to split a BAM per chromosome see http://code.google.com/p/jvarkit/wiki/SplitBam

It also creates an empty BAM (filled with a pair of mock SAMRecords) for each chromosome in the Reference, if no SAMRecord was found for the chromosome.
Question: Hello,

I did some googling about mRNAseq pipeline and built that pipeline. So I would like to know if my pipeline is okay.

Reference genome : Human gencode GRCh37
Aligner : STAR
Expression quantification : RSEM
Gene count table convert? : tximport
Statistical test : DESeq2

it seems like people use htseq or featureCount for gene count table and put it into edgeR or DESeq2. I am not sure, if I go for RSEM to run DESeq2 instead of htseq or featureCount. I tried kallisto and slueth, but I want to stick to this traditional alignment method.

Any advice would be apprieciated.
Thank you,

Answer: > I would like to know if my pipeline is okay.

RSEM can use STAR within its own script, so I would do that instead of running STAR yourself, unless you want to sophisticate your mapping more.

> I am not sure, if I go for RSEM to run DESeq2 instead of htseq or featureCount.

The output is comparable, it's a matter of choice, perhaps RSEM is more sophisticated. But in the end it's counting mapped reads per gene/transcript.

> I tried kallisto and slueth

Read this: http://www.rna-seqblog.com/limitation-of-alignment-free-tools-in-total-rna-seq-quantification/
Question: I have three samples each with 3 biological replicates i.e. one control and other at dose 12ug/ml and 25ug/ml at 24hr

```
> library('DESeq2')
> directory<-"./data"
> sampleFiles <- grep("Htseq",list.files(directory),value=TRUE)
> sampleCondition<-c("treated","treated","treated","treated","treated","treated","untreated","untreated","untreated")
> sampleGroup <-c ("A","A","A","B","B","B","C","C","C")
> sampleTable<-data.frame(sampleName=sampleFiles, fileName=sampleFiles, condition=sampleCondition, group=sampleGroup)
> sampleTable

sampleName fileName condition group
1 24h_12ug_1.count 24h_12ug_1.count treated A
2 24h_12ug_2.count 24h_12ug_2.count treated A
3 24h_12ug_3.count 24h_12ug_3.count treated A
4 24h_25ug_1.count 24h_25ug_1.count treated B
5 24h_25ug_2.count 24h_25ug_2.count treated B
6 24h_25ug_3.count 24h_25ug_3.count treated B
7 24h_ctrl_1.count 24h_ctrl_1.count untreated C
8 24h_ctrl_2.count 24h_ctrl_2.count untreated C
9 24h_ctrl_3.count 24h_ctrl_3.count untreated C
```

What should be good design if I want to see the differential expressed genes between treated vs untreated as well between different groups?

ddsHTSeq <- DESeqDataSetFromHTSeqCount(sampleTable=sampleTable, directory=directory, design=~group)

**OR**

ddsHTSeq<-DESeqDataSetFromHTSeqCount(sampleTable=sampleTable, directory=directory, design=~condition)

What will be the difference between the two designs?

I tried the following design but it gave me error

```
ddsHTSeq<-DESeqDataSetFromHTSeqCount(sampleTable=sampleTable, directory=directory, design=~group + condition + condition:group)
Error in DESeqDataSet(se, design = design, ignoreRank) : the model matrix is not full rank, so the model cannot be fit as specified. one or more variables or interaction terms in the design formula are linear combinations of the others and must be removed
```

How will the design will change if I add three samples each with 3 biological replicates i.e. one control and other at dose 12ug/ml and 25ug/ml at 6hr

Answer: Did you post this on SEQanswers too? I feel like I've seen it before.

If you want to look at treated vs. untreated and between group comparisons then just use ~group and use the appropriate contrast for the treated vs. untreated comparison.

The difference between ~group and ~condition is that in one case you're asking, "What's the difference due to treatment type?" while in the second case "What's the difference due to treating at all...ignoring that there are multiple treatments being used?" In reality, using ~condition will probably yield crappy results, since you have a batch effect that you're not compensating for.

You got that error because you can't simultaneously determine the effect of group C and the untreated condition...since they're the same. Similarly, the interaction terms are identical to the group terms, so they can't be determined either.

If you add more samples/conditions then the model will need to be changed according to the question(s) you want to ask. Given that you wanted inter-group comparisons before, I suspect you'll need to just add things as new groups and use appropriate contrasts for each question.
Question: ![enter image description here][1]

Check with this fig, I'm so confused about these

1. I think signal generation step is not belong to macs2, since it is
refer to mapping result, mapping let us know generated signal.

1. And peak calling is a function of macs2 but not mark here, actually,
I can't tell the difference between peak calling and signal
generation? For this fig, they need same input, bam or sam, right?

1. Which step's result can be as input to IGV? I used **fold change
over control** bigwig files as input, worked, but some slices says after peak calling?
2. **I just have fold change over control and signal p-value result
bigwig files now, so if I need to do peak calling, I must get
mapping results bam or sam? Can't transfer bigwig files as input, right？ I
mean information equivalence is not problem? Some tips, thanks！**

*I got figure from this website[https://www.encodeproject.org/chip-seq/transcription_factor/][2]*

[1]: https://gmgitx.github.io/blog/ENCPL493SGC.png
[2]: https://www.encodeproject.org/chip-seq/transcription_factor/

Answer: 1. You're right that the "signal generation" part is confusing there. MACS2 can output an optional bedGraph track (the `--bdg` option I think) that is essentially a coverage/signal track. It's highly likely that that's what's intended there. Truth be told I find bigWigs from `bamCoverage` to be more useful, but I'm rather biased.
2. Peak calling and signal generation are part of the same step, it's unfortunate that they show them to be unrelated, since you'd just run MACS2 once.
3. You can use the bigBed files and the bedGraph (N.B., you would be better served with a bigWig file, either by converting the bedGraph to bigWig or simply using `bamCoverage` from deepTools).
4. While it's possible to do peak calling on a bedGraph file with MACS2, I don't generally recommend doing it unless you really really know what you're doing. MACS2 `callpeak` essentially makes a bedGraph file internally, but you'd need to match the normalization it's doing and then run `macs2 bdgpeakcall`. I strongly encourage you to try and track down a BAM file instead.