title: Statistical test design question

text:
# Question
Title: Statistical test design question
I have condition as such mouse liver and spleen data. So in liver I have WT (Wild Type) and KO (knock out) as well as in spleen WT (Wild Type) and KO (knock out).

I have to compare the parasite load in WT and KO condition in both liver and spleen.

In terms of data points what I got is this:

* 10 mice for WT, 6 mice for KO in case of liver
* 10 mice for WT, 5 mice for KO in case of spleen.

```
WT 20 8 21 9 22 10 23 11 24
KO 12 25 26 27 28

```

So far what I know is I have to perform a parametric test.

But I'm not sure if I have to go for a paired t-test or an unpaired t-test.

Because here I'm comparing with between two groups, only so I don't have to go for anova if I'm not wrong. Secondly there is a difference in number of observation that is available for WT as compared to KO. So would that affect my choice of test as well.

Suggestion would be highly appreciated

# Answer
Here's a combination of general advice and my thoughts about your specific data:

If the samples are not matched, e.g. by time point, then it will be unpaired. You need a way to say sample 1 KO and sample 1 WT are uniquely matched to be paired, sample 2 KO and sample WT will be paired etc ... pairing is great for a time series. If you have replicates within the paired groups then ANOVA is the way to go, particularly if the sample sizes are uneven.

**Note**: in your case the sample sizes are uneven, you cannnot perform a paired T-test between two uneven sample sizes. HOWEVER you appear to a natural pairing, but I don't know whether this is biologically justified, i.e. liver and spleen. So you appear to have a biological pairing with uneven sample sizes. If this is correct it suggests a **1-way ANOVA** is a better idea than a T-test.

Parametric test: check for normality in each data set, e.g.

* A Q-Q plot against the normal distribution (a cool way to do this)
* Homogeneity of variances is Hartley's Fmax test (good for T-test and ESSENTIAL for 1-way ANOVA)
* 68% of the data should be less than ±1 standard deviation around the mean, 95% less than ±2 standard deviation and 99% less than ±3
* Kolmogorov-Smirnov and/or Shapiro-Wilk test
* Skewness & kurtosis magnitudes \< 1.96 respectively

**Note** if you combine liver and spleen data sets and find a nice normal distribution, and one or both is not normal by itself, then you would need to consider the biological justification for combining the data sets.

If the data passes all those tests, you are good to go. If not .... gets complicated,

* Normalisation, this is where you transform the data, in your case I'm not sure, but typically log transformation and then perform the checks for normality above. There are lots of ways to normalise data.
* Repeat the T-test/1-way ANOVA with normalised data
* Ignore it - yup its top level statistical advice, it is context specific and never a universally shared idea ... I suspect that is where the advice "just use a parametric test" came from. Please be aware if you do this, someone could turn around and sau "you've done it wrong". The statistical arguments about why this assumption can be valid get complicated.
* Use a non-parametric test, orthodox statistical advice, in this case it is the Mann-Whitney U test and **Kruskal–Wallis test** if you have paired replicates (i.e. non-parametric version of 1-way ANOVA).

Non-parametric tests often lack power and I suspect for this sample size that will be the case. Personally I would try and ignore all this advice :-) - okay I would definately check for normality and investigate avenues towards normality - then I would perform a general linear model (GLM). The reason is that it is a much more powerful method for identifying whether the two gropus are different. Performing a single GLM with two paired groups is not trivial, so I would have to give serious consideration to 1-way ANOVA, or KW-test, which would be a cleaner solution - but could miss a critical part of the data.
> 2 votes

# Answer
> So far what i know is i have to perform a parametric test .

How did you determine this?

> But im not sure if I have to go for a
>
> Paired t-test or Unpaired t-test

If you can't tell the difference between when to do a paired t-test and an unpaired one, you need to stop and talk to someone who can teach you basic statistics. Throwing statistical tests at numbers when you don't understand them is a recipe for disaster.
> 2 votes

---
Tags: statistics
---


title: rentrez to retrieve NIH files

text:
# Question
Title: rentrez to retrieve NIH files
I'm trying to use rentrez package to retrieve NIH files, and ran the following code (I followed what was on the rentrez tutorial):

```
pubmed_search <- entrez_search(db = "pubmed",
                    term = "Case Reports[Filter] AND cardiovascular disease AND English[lang] AND 2009:2019[PDat])",
                    retmax = 792711,
                    use_history = TRUE)

for(i in seq(1,10000,100)){
    recs <- entrez_summary(db="pubmed",
                           web_history=pubmed_search$web_history,
                           retmax=100,
                           retstart=i)   cat(
seq_start+99, "sequences downloaded\r") }

```

However, instead of the 10,000 texts I expected, I only downloaded 100 files. Could so

# Answer
> 1 votes
I don't program in R, but replace the last line of you code with the following and I suspect you will find it works,

```
for(i in seq(1,10000,1)){ recs <- entrez_summary(db="pubmed", web_history=pubmed_search$web_history, retmax=100, retstart=i) cat(seq_start+99, "sequences downloaded\r") }

```

**Rationale** I would guess `seq` is similar to the `linspace` function in Python, and put simply `10000 (total sample size) / 100 (the step size) = 100`, 100 being the number of downloads you observed. In the code above its `10000/1`. I would also suspect "step size = 1" is default, if this was correct then `seq(1,10000)` would give you your required output.

As a random aside I do like the presence of curly brackets defining the loop, it reminds me of Perl.

---
Tags: r, entrez
---


title: How to create a nexus file for BEAST?

text:
# Question
Title: How to create a nexus file for BEAST?
I am learning to use BEAUTI and BEAST software. I followed their tutorial and was able to run example data without any problem. Now, I want to do phylodynamic analysis on my data. I have downloaded 70 picornavirus sequences from the gene bank. My questions are :

1. What tool do I have to use to create Nexus file that is used by BEAUTI ?
2. How do I enter the information about the date (year ) the virus was isolated and location information which is required for rates and dates (evolutionary) analysis in BEAST, spatial analysis by SpreaD3 ?

---

Notes on SpreaD3 is a update on SPREAD by Bielejec et al. SPREAD: Spatial Phylogenetic Reconstruction of Evolutionary Dynamics. Bioinformatics, 2011, 27:2910-2912. This uses spatial Bayesian reconstruction based on the phylogenetic output of Beast

# Answer
Its easy,

1. Download the data from NCBI in fasta format (its an option in one of the pull down headers).
2. Load it into the latest version of MEGA, version 7(?), perform the alignment therein.
3. Export from MEGA as Nexus format and upload it into Beast

You could of course perform the alignment in Clustalx, Clustalo or Muscle or T-Coffee, (hmmm I could write an alignment program just to name it T-4two (sorry British joke)). If the output option of "nexus" was not present in the alignment program of your choice, then again upload the alignment into MEGA and export it as nexus.

Anyway, I could go into specific detail about the various points and clicks therein, but I think its easy to figure out. Most phylogeny programs will do an automatic fasta-\>nexus conversion, but anyway.

Its easy, Beast is not trivial to use properly however and ensuring convergence is not at all easy, in addition the multiple simulation a priori models, whilst powerful, give different answers. To try and explain this a different way, it is easy to get an answer in Beast, but difficult to get a robust answer from it and thats the bit I would give serious consideration about. It does amaze how widespread the program is (I guess its the nice GUI), given how complicated it is to obtain robust result.

---

Just to add you could try the online tool here for the fasta to nexus conversion: . I've not tried it, its just an after thought. If you do Python, Biopython will do it, but its just as easy to code.
> 1 votes

---
Tags: sequence-alignment, ngs, covid-19, virus
---


# Question
Title: How can I translate infix to mathml using the python interface to libsbml?
Assume I have an infix expression such as k1\*s1, how do I use libsbml to translate this expression to mathml using the python interface?

# Answer
It’s quite straightforward to convert infix to mathml using libsbml. Make sure you have libsbml installed in python by typing

```
pip install python-libsbml

```

at the windows/Mac/Linux command line or

```
!pip install python-libsbml

```

if you are running a Jupyter notebook. Some tools such as Tellurium already have libsbml preinstalled.

At the python command line type

```
import libsbml

```

to import the package. To enter the infix expression type the following:

```
p = libsbml.parseFormula ('k1*s1')

```

To generate the mathml type

```
print (libsbml.writeMathMLToString (p))

```
> 4 votes

---
Tags: modelling, systems-biology, sbml
---


# Question
Title: Looking for a taxonomic dataset that links tree genera (or species) to families
I am looking for a taxonomic dataset that would allow me to links tree genera (or species) to families. In particular, I am interested in the tropical species found in the Amazon.

I found a long but not exhaustive list on Wikipedia: https://en.wikipedia.org/wiki/List\_of\_trees\_and\_shrubs\_by\_taxonomic\_family

After processing it, I realised it contains information about 404 genera, while the dataset I am dealing with has around 600. (in the end, I only got 23% of my data covered, very sad)

Does anyone know if such a crosswalk file exists, in some scientific article or somewhere in the Internet?

Moreover, please tell me if there are other StackExchange websites which could be more appropriate, I wasn't sure of where to post this.

# Answer
Probably the best currently available online database is RBG Kew's Plants of the World Online. It organises taxonomic, morphological and distribution data from a range of databases. You can search for your taxa of interest and find their associated higher ranks.

Alternatively, there is World Flora Online, replacing the now-defunct The Plant List v1.1, which also integrates taxonomic data from a range of reliable sources and national lists.
> 2 votes

# Answer
The NCBI taxonomy is a quite extensive resource: https://www.ncbi.nlm.nih.gov/taxonomy

I would thus expect there is a good chance that it covers all of the taxa that you are interested in. It officially claims that it is not an authoritative source, but it's still pretty darned good.
> 4 votes

---
Tags: database, public-databases, taxonomy
---


# Question
Title: Is there a way to use BioLemmatizer from Python?
Lemmatization is one of the most important tasks for text normalization in natural language processing (NLP), i.e. using deep learning to understand text data.

One of the good lemmatizers I found for biomedical domain is BioLemmatizer. However, it is developed in Java and is available from maven repo. I would like to use it from Python, but there seems to be no counterpart for BioLemmatizer in Python. Is there a way I can access BioLemmatizer from Python using some module like JPype. JPype, however, supports access only for the default classes in Java.

# Answer
This area is huge in Python because of machine learning text processing, in fact O'Reilly have an entire book on it . However, we don't really deal with words in Bioinfo SE,

You might want to look at the nltk (natural language toolkit) for Python, but for direct parallels with the Java package bio package (your question) I don't think it's something we would know.

---

You could solve this via subprocessing into Java, and passing the result back to Python. I get the idea that the core ML/ DL used in natural language is set up for Python. Google for example have fantastic tensorflow stack for this.
> 0 votes

---
Tags: python, text-processing, java
---


# Question
Title: Filter reads belonging to unique sequences with threshold
I have CLL samples as fastq files and I want to remove those reads which have a unique sequence with less than 10 read counts each. I tried this by following way using awk to make it faster:

```
for fq in $(ls fastq1/*.fq.gz); do
    basename=$(basename "$fq" | cut -f1 -d '.')
    zcat $fq | awk "/^@/{getline; print}" |sort | uniq --count | awk '{if($1>= 10) print $2}' > tmp.txt; grep -B 1 -A 2 -f tmp.txt $fq | grep -v ^--'| pigz -p 10 > fastq2/$basename.fastq.gz

done

rm tmp.txt

```

But this is not working well as many reads with more than 10 read counts are also missing. I think it's OK to stick with bash/awk for fast fastq file processing. What am I missing here, or is this command wrong. Can anyone suggest me the right and fast way to do it??

# Answer
Your algorithm generally works, you just had one minor issue, and a few stylistic changes that could be made.

I'd adjust your script to the following:

```
for fq in fastq1/*.fq.gz; do
    basename=$(basename "$fq" .fq.gz)
    zcat "$fq" | awk "/^@/{getline; print}" | sort | uniq --count | awk '{if($1>= 10) print $2}' > tmp.txt
 zcat | grep --no-group-separator -B 1 -A 2 -f tmp.txt "$fq" | pigz -p 10 > "fastq2/$basename.fastq.gz"
done

rm tmp.txt

```
1. Instead of

`for fq in $(ls fastq1/*.fq.gz)`

use this

`for fq in fastq1/*.fastq.gz`

Iterating over `ls` is generally not ideal, and can give unexpected results. Better to just rely on bash globbing.

2. Purely stylistic, but I'd simplify and use `basename` to remove the suffix, instead of `cut`.
3. Rather than removing the separators that `grep` will introduce by filtering them out, use the following option to not create them: `grep --no-group-separator`
4. And, minor bug, as others have mentioned, you should use `zgrep` or `zcat | grep` since your files are compressed.

Oh - and quote your variables.

Finally, while this does work, it may not be the most efficient method. The file is being decompressed twice, and may be better to decompress once to a temporary file that can be used and then deleted.
> 1 votes

# Answer
Checking for line start to find the read header is not the best option. Even if unlikely it is possible that the qualitys starts with a @. It's better to compare the current line number. Also you can use awk directly to count and print the result.

```
zcat $fq|awk 'NR%4==2 {count[$0]++} END { for(sequence in count) { if(count[sequence]>=10) print sequence }}'

```

As the other statet out, the second command should be fine if you use `zcat` instead of `cat`.
> 1 votes

---
Tags: fastq, bash, awk
---


# Question
Title: How to remove degenerate sequence (N4) at the 5' and 3' of the reads?
I have to analyse GSE114327. For that, I use GSE instruction based on GSM3139597.

In that data processing section has written:

> FASTQ reads were subjected to adapter removal (cutadapt) per manufacturer's instructions (NEXTflex Small RNA Sequencing Kit v3).

So, I use `cutadapt` for trimming. Then

> The degenerate sequences (N4) at the 5' and 3' end of the reads were removed

based on the instruction. But I don’t understand that sentence and I don’t know how I can remove degenerate sequence (N4) at the 5' and 3' of the reads.

# Answer
From the link you posted it looks as though they used UMI tools as per the section which says *'were removed and added to the name of the read as a unique molecular identifier (UMI) using UMI tools'*. The link above has instructions for installation and use. From this I would guess that this data has Unique Molecular Identifiers (UMIs) for each read.
Looking up the protocol for the nextflex kit used (See Link) confirms this:

> **Data Analysis** The 3' and 5' adapters included in this kit both contain 4 random bases that will appear immediately 5' and 3' to the insert in sequencing data. The presence of these random bases should be considered when choosing an alignment strategy. When using “end-to-end” alignment, we recommend processing data in the following manner: 1. Clip the 3' adapter sequence (TGGAATTCTCGGGTGCCAAGG). 2. Trim the first and last 4 bases from the adapter-clipped reads. 3. Perform alignments as normal.

It seems that N4 is a 4 base random barcode at both the 5'& 3' end of the read - these are cut from the read and this metadata is added to the read name so that every read can be uniquely identified.

*Edit:*

You can install UMI tools via:

`conda install -c bioconda -c conda-forge umi_tools`

or

`pip install umi_tools`

Usage instructions can be found at:

https://umi-tools.readthedocs.io/en/latest/reference/extract.html

> Usage: For single ended reads, the following reads from stdin and outputs to stdout:
>
> umi\_tools extract --extract-method=string --bc-pattern=\[PATTERN\] -L extract.log \[OPTIONS\]
>
> For paired end reads, the following reads end one from stdin and end two from FASTQIN and outputs end one to stdin and end two to FASTQOUT:
>
> umi\_tools extract --extract-method=string --bc-pattern=\[PATTERN\] --bc-pattern2=\[PATTERN\] --read2-in=\[FASTQIN\] --read2-out=\[FASTQOUT\] -L extract.log \[OPTIONS\]
> 3 votes

---
Tags: fastq, data-preprocessing
---


# Question
Title: semantic similarity measurement for cell line ontologies
I have a set of cell line pairs and I want to know to what extent the pairs are similar based on their ontologies. The problem I have is that I have found a Python library called Fastsemsim, but it has many similarity measurements. Or to be exact, there are actually many ontology similarity measurements. And I want to know which similarity measurement is best suited for my case (cell line ontologies similarities). There might be a couple of good ones that are good for cells ontologies. These are a set of similarity measurements that are available in Fastsemsim.

Semantic Similarity

* SemSim.TermSemSim
* SemSim.ObjSemSim
* SemSim.ObjSetSemSim
* SemSim.SetSemSim
* Specific Semantic Similarity measures
* SemSim.MixSemSim
* SemSim.avgSemSim
* SemSim.maxSemSim
* SemSim.BMASemSim
* SemSim.SemSimUtils

Source: https://fastsemsim.readthedocs.io/en/latest/library.html

# Answer
> 1 votes
About the metric you use, there are several reviews on the topic, if I recall well I think that the information content based are the best ones, of those, I only recognized the JiangConrathSemSim listed on the webpage.

However, if you want to use a semantic similarity measure I would use the BMA method, which is more comprehensive than other methods. But the BMA method is only a way to aggregate several metrics of ontologies into a single numeric value. I am not confident that FastSemSim implements correctly those methods.

Note, that according to the license, you should contact the authors before using it:

> You must contact us before using it or its results or any application/tool based on it in any published work.

---

I am not sure that the cell lines fit in a ontolgogy. Because ontologies are defined as directed acyclic graph (DAG) and there are some cycles that happens in cells like the induced pluripotent stem cells show. However they might be right if you just compare the cell lines and not the state of the cells on those lines (use it only if you want to compare the lineage of your cells lines).

---
Tags: python, single-cell, cell-line, ontology
---


# Question
Title: Hierarchical clustering for outlier detection - single cell RNASeq & WGCNA
I wish to conduct WGCNA on a single cell RNASeq dataset and, when choosing the optimal beta parameter after running `pickSoftThreshold` am presented with a rather strange plot:

I have normalized the data using `NormalizeData` in `Seurat` using the LogNormalize method.

One of my colleagues suggested using `hclust` as a method for outlier detection to see potential cells/samples that could be excluded. I think this makes sense however I haven't found much information on such procedures besides the following paper: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5845381/.

My question would be, is this procedure just be as simple as running a hierarchical clustering on the data and seeing which samples cluster by themselves (and then removing those based on a minimum cluster size criteria) or is their more to it? Also, for those familiar with WGCNA, why does the scale-free topolgy plot generate such strange looking values, is this indicative of outliers being present in the data?

Thank you for your time.

# Answer
For outlier identification I suggest using the sample network approach developed by Oldham et al. It basically amounts to constructing a inter-sample connectivity (assuming `normalizedData` contains the log-normalized data with genes in colunms and samples in rows)

`k = colSums(cor(t(normalizedData)))`,

scaling the connectivities

`Z.k = scale(k)`

and then removing samples with high negative Z.k (e.g. Z.k \< -4).

A slightly different approach is based on euclidean distances rather than correlations,

`Z.k = scale(as.matrix(dist(normalizedData)))`

and removing samples with high *positive* Z.k (since distance is a measure of dissimilarity).

Regarding the scale-free topology fit plot, I wouldn't worry about it too much, the fit index sometimes jumps around a bit. It certainly does not look like you have strong outliers or very heterogeneous data. But it is always a good idea to plot the sample clustering tree to check (1) that there are no strong clusters, and (2) that the outlier identification did something reasonable.
> 2 votes

---
Tags: r, scrnaseq, seurat, clustering, wgcna
---


title: Hard-clip primers in a bam file

text:
# Question
Title: Hard-clip primers in a bam file
I am working on a sample which was sequenced using an amplicon kit and I would like to ensure that no variants are called in the PCR primer sequences. I know that the primer sequences are part of the genome and a variant there is no less valid than a variant anywhere else, but for reasons I cannot control, I need to ensure that no variants can be called in the primer regions.

My first thought was to use a soft-clipping tool like BAMClipper. However, the variant caller I am using (Sentieon's implementation of GATK) preforms a local realignment when calling and if it finds that soft-clipped regions align perfectly against the genome, it will remove the soft clips. This makes perfect sense for cases where the softclips were added by the aligner, but sadly breaks my workflow since I need the softclipped primers to remain clipped.

I would therefore like to hard clip, to completely remove the primer sequences from the bam file before variant calling. I have found one tool that seems to be able to do this, TrimPrimers by Fulcrum Genomics, but this fails on my input bam file with an error I can't really make much of.

Is there any other tool that can take a list of regions (e.g. a bedpe file) and a BAM and hard-clip the regions from the BAM?

# Answer
Another option would be to adjust the quality of the primer bases to be below the threshold that your variant caller uses. In my case I have "masked" the primer bases by dropping the quality to 0. This prevents the variants from being called based that are synthetic but if you have an overlapping amplicon it will call bases from those reads.
> 2 votes

# Answer
> I know that the primer sequences are part of the genome and a variant there is no less valid than a variant anywhere else,

You won't see variants under PCR primers. You get the primer sequence itself in the vast majority of reads.

I wouldn't clip the bam file, I'd clip either the fastq, or clip within bcl2fastq while making the fastqs.
> 1 votes

---
Tags: sam, primer, base-clipping
---


# Question
Title: Compare functional enrichment between genelists
I want to compare the results from functional enrichment (e.g GO term enrichment generated by gProfileR) between two gene lists. The idea is to see if these gene lists are enriched in the same terms or if they have specific terms.

This is easy to do at a glance when you have few terms but when the number of terms increases and due to redundancies in GO term annotation a manual comparison becomes unreliable.

So the questions is: **Is there a tool for comparing the results from functional enrichment between different gene lists?**

# Answer
> 1 votes
To compare gene list based on gene ontologies you can use GOSemSim (Note, I'm a contributor to the package). You can select the semantic similarity to compare them, that is how the redundancies and hierarchical relations are counted. The main result you can obtain is similarity matrix for each comparison between gene lists.

However if you perform a GO term enrichment and then you compare the gene list of each GO term by looking at which GO do they belong it won't be much useful. I would suggest you to use BioCor (I'm the developer of it), which can use any pathway information (such as the ones of Reactome) to asses if they are related to the same pathways.

---
Tags: functional-annotation, go-enrichment
---


# Question
Title: For phylogenetic tree construction from core-genome which one is preferable: amino-acid based MSA or nucleotide based MSA?
The genomes are from same species. Is it true that, in phylogenetic tree constructed from amino-acid based MSA (multiple sequence alignment) some information are lost, so for phylogenetic reconstruction of closely related sequences nucleotide based MSA should be used?

# Answer
For a species complex nucleotide phylogenies all the way.

The reason is neutral mutation, which is not observed at the amino acid level because these mutations are part of the degenerate code which don't induce amino acid substitutions. Such mutations are commonly referred to as the third codon wobble, although occasionally the first codon can mutate and not produce a change in amino acid.

Neutral otherwise known as synonymous mutations occur at a much higher rate than mutations that cause a substitution of an amino acid. This is the reason why you were told information is lost at the amino acid level because you're losing neutral mutations.

There are circumstances when the amino acid phylogenies are preferred to nucleoside trees, but these are bit pretty specialised areas.
> 4 votes

---
Tags: sequence-alignment, fasta, phylogenetics, phylogeny
---


# Question
Title: Paired end sequencing with R1 and R2 of different length. Possible?
There are commercial sequencing kits/sequencers that allow for paired end sequencing in which the two reads obtained for each fragment are of different length?

For example 150bp for R1 and 50bp for R2?

Is it at least possible in principle? There are some published speculations about that possibility?

# Answer
As indicated in the comments, yes, you can definitely tell standard Illumina sequencers to sequence mates in a pair to different lengths. This is quite common in single-cell RNA-seq, where R1 is often just a cell barcode and an UMI, so sequencing it to 150 bases would be a huge waste or reagent. One caveat to this is that you can't buy reagents sized for these custom length, just things like 150 base PE or other standards. So unless you end up saving and pooling reagents after such runs you end up paying the same for projects with one shorter read than if you'd just directly sequenced everything to the same length. Given that sequencer runs only take a day or two at the moment, this also doesn't result in that immense of a time savings.
> 5 votes

---
Tags: ngs
---


title: help with CRCmapper

text:
# Question
Title: help with CRCmapper
I'm using CRCmapper on my data. There are several modules in the tool. I've used exactly instructed in the usage.

Commands used:

```
python CRCmapper.py -e MB12_noEncBlack_AllEnhancers.table.txt -b MB12.bam -g HG19 -f hg19-stdChr/ -s MB12_peaks.bed -n MB12 -o MB12/

```

Here is the error I got:

```
Warning: text mode turns off computation of q-values
IDENTIFY TF-TF INTERACTIONS
Traceback (most recent call last):
  File "CRCmapper.py", line 605, in <module>
    main()
  File "CRCmapper.py", line 593, in main
    graph = buildNetwork(projectFolder, projectName, candidateGenes, refseqToNameDict, motifConvertFile)
  File "CRCmapper.py", line 381, in buildNetwork
    target = refseqToNameDict[region[0]]
KeyError: 'AP1'

```

My expertise in python is very basic. I would be grateful if someone take a look at this error and suggest how to fix.

P.S. I contacted authors a long back, I didn't get any response.

Thank you.

# Answer
It seems that the author of CRCmapper.py wants to extract `sequence_name`, `start`, `stop` columns from `fimo.txt` file. But the code extract `motif_alt_id`, `sequence_name`, `start` columns. So the index of the code should be changed.

`fimo.txt` looks like this:

```
motif_id    motif_alt_id    sequence_name   start   stop    strand  score   p-value q-value matched_sequence
Transfac.V$AP1_02   AP1 NM_003046931|3|95762|152118 2332    3118    -   11.9802 6.49e-05        TGACTCA

```

In line 377-382 of the code:

`line[1]` should be changed to `line[2]`.

```
for line in fimoTable[1:]:
    source = motifDatabaseDict[line[0]]
    region = line[1].split('|')
    target = refseqToNameDict[region[0]]
    location = (region[1], int(region[2]), int(region[3]))

```

In line 391-394 of the code:

`line[2]` should be changed to `line[3]`, `line[3]` should be changed to `line[4]`.

```
    # Count unique motifs
    if (region[1], int(region[2]) + int(line[2]), int(region[2]) + int(line[3])) not in motifDictSE[source]:
        edgeCountDictSE[(source, target)] += 1
        motifDictSE[source].append((region[1], int(region[2]) + int(line[2]), int(region[2]) + int(line[3])))

```
> 2 votes

# Answer
This error came from the fact that the output of the new fimo version has an additional 'motif\_alt\_id' column in second position.

CRCmapper has been modified to take this into account. You can find the updated version here: https://bitbucket.org/young\_computation/crcmapper/src/master/

P.S.: sorry if we missed you request user1545 and thanks for using CRCmapper!
> 2 votes

---
Tags: python, networks
---


# Question
Title: In SBOL, when should I annotate a DNA sequence vs. making sub-components?
When representing a DNA sequence in SBOL as a `ComponentDefinition`, you can mark things like promoters and coding sequences in two different ways:

* either as a `SequenceAnnotation` in the `ComponentDefinition`, or
* as a `Component` linking to a `ComponentDefinition` for the sub-comoponent

When should you choose one over the other?

# Answer
Typically you want to use a `ComponentDefinition` when an sub-structure is something you might want to pull out and re-use in another genetic design, and a `SequenceAnnotation` whenever it is not.

For example:

* Promoters, terminators, and protein coding sequences are typically best represented with their own `ComponentDefinition`, since these are often pulled out and re-combined in new genetic designs. Representing such an element with a `ComponentDefinition` makes it easy to find the commonalities between designs.
* Assembly scars, binding sites, and cut sites are typically best represented with a `SequenceAnnotation` because they are typically not very interesting except in the context of their surrounding genetic context. Representing such an element with a `SequenceAnnotation` helps keep them associated tightly with that context.

Not every case is cut and dried, however. For example, ribosome entry sites are often treated as separable components, but their performance is very tightly tied to the specifics of the coding sequence that they modulate. Thus, in some cases it may make more sense to represent them separately as a `ComponentDefinition` and in others to represent them as a `SequenceAnnotation` on a joint RES/CDS `ComponentDefinition`.
> 4 votes

---
Tags: sequence-annotation, sbol
---


# Question
Title: How do I select a proper design matrix to use in DRIMSeq?
I am trying to use DRIMseq for DTU with 2 treatments on two different strains of an animal model. How do I select a proper design matrix to use in DRIMSeq ? How I know that one is the correct one to use? And how do I see which transcripts are differentially used among conditions after the pipeline was run? I suspect it is in the plotProportions function.

Also, where do I see the logFoldChange between conditions?

Many thanks in advance as it would help me a lot in my work.

```
> sessionInfo()
R version 3.5.3 (2019-03-11)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] sva_3.30.1          BiocParallel_1.16.6 genefilter_1.64.0   mgcv_1.8-28
[5] nlme_3.1-137        edgeR_3.24.3        limma_3.38.3

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.1           compiler_3.5.3       pillar_1.4.2         bitops_1.0-6
 [5] tools_3.5.3          digest_0.6.20        zeallot_0.1.0        bit_1.1-14
 [9] annotate_1.60.1      RSQLite_2.1.1        memoise_1.1.0        tibble_2.1.3
[13] lattice_0.20-38      pkgconfig_2.0.2      rlang_0.4.0          Matrix_1.2-17
[17] DBI_1.0.0            rstudioapi_0.10      yaml_2.2.0           parallel_3.5.3
[21] IRanges_2.16.0       S4Vectors_0.20.1     vctrs_0.2.0          locfit_1.5-9.1
[25] stats4_3.5.3         bit64_0.9-7          grid_3.5.3           Biobase_2.42.0
[29] AnnotationDbi_1.44.0 survival_2.44-1.1    XML_3.98-1.20        blob_1.2.0
[33] matrixStats_0.54.0   backports_1.1.4      splines_3.5.3        BiocGenerics_0.28.0
[37] xtable_1.8-4         RCurl_1.95-4.12      crayon_1.3.4

```

# Answer
> 2 votes
From your question I understand that you have two treatments in two mouse strains.
If you are interested in general treatment effects independet of the animal strain you would use the both variables in the design formula (similar to controlling for a batch effect).

```
meta.data <- data.frame(treatment = c("a", "a", "b", "b"), strain = c("x", "y", "x", "y"))
design_combined <- model.matrix(~ treatment + strain, data = meta.data)

```

If you want to make the comparisons for each strain separately, the easiest way will be to split your data set by strain and do a set up just for treatment. Alternatively, you can also merge the two columns into a new one in case you also want to also compare between the strains.

```
meta.data$combined <- paste0(meta.data$treatment, "_", meta.data$strain)
# refactor to change the reference level
meta.data$combined <- factor(meta.data$combined, levels = "bx", "by", "ax", "ay")
design_combined <- model.matrix(~ combined, data = meta.data)

```

To get the fold changes, you can extract the proportions for each transcript using

```
DRIMSeq::proportions(dataSet)

```

And from the proportions the log2 fold changes can be calculated with ease.

However, when I look at your plot something seems very off with your analysis right now. In all the examples I have seen and also in my own data the diamonds are representing the mean proportion of a transcript in a condition. See this example from the DRIMSeq manual on Bioconductor:

In your case there are several different means per group and if you look at the line to the very right the two treatments seem to be to be right on top of each other.
So I think that treatment and strain are somehow mixed in your setup and this is causing this kind of merged plot. If you use a design as outlined above you should get much better results. Additionally, you could maybe check that the meta data is correct.

Also in case you have not read it yet, the manual is a great resource:
https://bioconductor.org/packages/release/bioc/vignettes/DRIMSeq/inst/doc/DRIMSeq.pdf

---
Tags: r, rna-seq, design
---


# Question
Title: Biological replicates on Chip-seq Transcription factor data
I have 2 biological replicates for chip-seq transcription factor data and I have applied macs2 peak calling for each replicate separately.

How can i make use of the biological replicates to extract the most out the data?

Thanks!

# Answer
You could first look at the degree of correlation between the two replicates - what proportion of peaks are shared between the two, versus peaks found in one sample only? This will give an idea of how repeatable the analysis is, and how many peaks are variable between samples or due to artefacts of the method.

When it comes to interpreting the biology behind the binding patterns, you are probably going to want to only interpret peaks which you find in both replicates.
> 2 votes

# Answer
This question is somewhat generic, so a generic answer is that ENCODE has a Transcription Factor ChIP-seq Data Standards and Processing page that can give you a useful starting point.

For TF ChIP-seq data with replicates, the Irreproducible Discovery Rate (IDR) method helps leverage replicates to produce higher confidence peak calls, producing both "optimal" and "conservative" peaks, defined by some IDR threshold. Here's ENCODE's description of the method:

> A statistical procedure that operates on the replicated peak set and compares consistency of ranks of these peaks in individual replicate/pseudoreplicate peak sets. Peaks with high rank consistency are retained. IDR can operate on peaks across a pair of true replicates resulting in a “conservative” output peak set, or across a pair of pseudoreplicates resulting in an “optimal“ output peak set. Peaks in the conservative peak set can be interpreted as high confidence peaks, representing reproducible events across true biological replicates and accounting for true biological and technical noise. Peaks in the optimal set can be interpreted as high-confidence peaks, representing reproducible events and accounting for read sampling noise. The optimal set is more sensitive, especially when one of the replicates has lower data quality than the other.

The Python code to run IDR on peak sets can be found on ENCODE's GitHub Chip-seq pipeline repo.
> 2 votes

---
Tags: chip-seq, macs2, transcription-factors
---


# Question
Title: making a customer reference by using cell ranger
Anyone here knows how to add a new gene (GFP) to customer reference by using cellranger? We use the Cre/loxp system by using R26RmTmG mice and want to use GFP to isolate specific cells. I have read the cellranger workflow but still, have no idea how to do it. Like it is said that I should update the GTF file as shown in the screenshot.

# Answer
> 1 votes
Step 1: Add the the GFP sequence to your fasta file

```
>GFP_name
ATG.....TAA

```

Step 2: Annotate the gtf so that it contains all the necessary information to annotate your information

Add this line to your GTF file but modify where necessary: "Make sure that all columns are tab separated"

```
GFP_name    AddedGene   exon    1   1966    .   +   1   gene_id "GFP"; transcript_id "GFP"; gene_name "GFP"; transcript_name "GFP"

```
* The first column has to correspond to the header of the fasta sequence
* start and end indicates the length of your sequence. I added a 1966 bp long sequence to my gtf

Code example: step 1:

```
cat reference.fa gfp_seq.fa > new_reference.fa

```

step 2:

```
cat genes.gtf gfp.gtf > new_genes.gtf

```

\- gfp\_seq.fa is your custom gene in fasta file \- gfp.gtf is your custom gtf line or lines if you want to add more

You can also open your current reference fasta and gtf file in a text editor and copy paste the sequence and gtf line in there.

Then run `cellranger mkref` on your new\_reference and new\_genes file

# Answer
> 0 votes
Update your gtf file, update the genome fasta, run cellranger mkref.

---
Tags: single-cell
---


# Question
Title: Does SBOL require representing intermediate products like mRNA?
Do I need to include intermediate gene products like mRNA in order to express a gene regulatory relationship in SBOL?

For example, let's say that I'm representing transcriptional regulation of the pTet promoter by means of a coding sequence that expresses TetR. Do I need to represent the TetR transcript mRNA as well as the TetR protein?

# Answer
> 2 votes
As far as I understand, the purpose of SBOL is not to create a comprehensive model, but to *communicate a design*. So the question you have to ask is "is the mRNA part of my design, or is it an implicit part of the biology which does not help to convey what I'm trying to build?".

I think your summary is a giveaway:

> For example, let's say that I'm representing **transcriptional regulation** of the **pTet promoter** by means of a **coding sequence that expresses TetR**.

Notice that the mRNA was not mentioned? I also daresay that if you drew this design as a diagram, the mRNA would be omitted - another indication that it is not helpful in communicating what you designed.

---
Tags: modelling, sbol
---


# Question
Title: Why does SBOL have both Component and ComponentDefinition?
What is the point of having both a Component and a ComponentDefinition in SBOL? Shouldn't you be able to link straight from one component to its sub-components?

The same question goes for Module and ModuleDefinition.

# Answer
> 2 votes
The key to understanding this is to think about a system with more than one copy of a genetic construct.

For example, let's say you have a construct that uses two copies of the pBAD promoter, each regulating a different coding sequence. Both of these promoters are pBAD promoters, with the same sequence and type information, so they should share the same `ComponentDefinition`.

But there are two copies. In the context of the larger construct, we need to be able to talk about each copy separately. For example, we might need to say that the first comes after a MoClo cut site, while the second comes after a terminator. For that, we need the `Component` object, to have a way to say "this particular copy of pBAD", as opposed to "pBAD in general", which is what the `ComponentDefinition` tells us about.

# Answer
> 1 votes
Essentially, a `ComponentDefinition` represents *the component itself*, while a `Component` represents *the component being used as a sub-component*.

There is a proposal within the SBOL community as to whether for SBOL version 3 `ComponentDefinition` should be renamed to `Component` and `Component` should be renamed to `SubComponent` to make this situation less confusing.

---
Tags: sbol
---


title: Extract mapping coverage from GTF files

text:
# Question
Title: Extract mapping coverage from GTF files
I am doing paired end RNA-seq analysis. I used STAR mapper and then StringTie for getting quantification of my data. After successfully running StringTie, I got `.gtf` files which will give me information about genes, their coverage along with TPM and FPKM values. It also produced some more `.ctab` files.

I want to know how what percentage of my reads mapped to reference genome. I know this can be done with BAM or BED files produced after STAR alignment. However, I want to know if this can be done with output files generated by StringTie?

I just want to clarify my question after reading comments and answers I am getting. I want to know if this can be done with `.gtf` files of StringTie. As I mentioned before, I already know how to achieve this with BAM files. This is little bit weird situation because, I have deleted my BAM files to get more space on my hard drive. So I could just run entire analysis again, but wanted to know if it is possible without redoing whole analysis. Apologies for the confusion.

# Answer
So consider this, you know how many reads map to each of your genes (1). Does that tell you anything about (2) the number of reads that mapped against non-gene regions (3) the number of reads that did not map against the genome at all?

If you only use the gtf files, you know the number of (1) but not of (2) and (3). To me it sounds like you would want to do something as described in this post.

What it boils down to is:

Finding out how many (unique) reads are mapped, counting pairs once

```
samtools view -F 0x4 foo.sorted.bam | cut -f 1 | sort | uniq | wc -l

```

Counting how many total reads you have. Since you have paired-end data the number of forward reads suffice:

```
cat forward_reads.fastq | grep '@' | wc -l

```

The percentage is then simply `mapped_reads/total_reads`
> 1 votes

# Answer
> I want to know how what percentage of my reads mapped to reference genome.

And the Log.final.out file from STAR doesn't tell you that?
> 0 votes

---
Tags: coverage, star, stringtie
---


# Question
Title: Single Cell RNA seq Analysis for low input Cells?
I am having a single cell RNA seq data from ~500 Cells. I do understand this number is very low for present day analysis tools like Seurat. So, I want to know what if one have a such a small number of cells which tool works best?

Thanks

# Answer
You can perfectly use any of the single cell processing packages out there. Seurat, Monocle3, Scanpy, Kallisto, etc. Following some of the tutorials and vignettes will most likely help you figuring out how to proceed.

These packages provide useful functions to explore single cell data, process the data and each of them has a series of advanced features to address specific questions. You don't need 10k cells to look at number of gene/cell, transcript/cells, % mitochondrial genes, etc. You can perfectly normalize and standardize 500 cells and run PCA, and analyze the complexity that describes your data.

To address if 500 cells is enough for you pilot we need to know what question you want to address. The following tool for example gives some intuition on how many cells to sample if you wish to detect n-cells per cluster.

https://satijalab.org/howmanycells

Hopes this helps
> 3 votes

---
Tags: rna-seq, scrnaseq, single-cell
---


# Question
Title: How do I add external labels and connecting lines to a Ternary Plot in R?
I have plotted a triangular plot,that three dimensions of which represent three herbal strategies. One dimension represents the strategy of C (competitive plant), the second dimension “S” (stress tolerant plants) and the third dimension”R” (ruderal plants), as you can see in the image below, the points on it represent the plant species. Now I am trying to write name of the plant species outside the triangle and connect with the arrow to the corresponding point. How do I do this work? Many thanks in advance as it would help me a lot in my work.

```
TernaryPlot()
#Plot two stylised plots side by side, and plot data
par(mfrow=c(1, 1), mar=rep(0.3, 4))
TernaryPlot(atip='C%', btip='R%', ctip='S%',
            point='UP', lab.cex=0.8, grid.minor.lines=0,
            grid.lty='solid', col='#FFFFFF', grid.col='GREY',
            axis.col=rgb(0.1, 0.1, 0.1), ticks.col=rgb(0.1, 0.1, 0.1),
            padding=0.08)
data_points <- list("Bromus dantonia" = c(47,   59, 149),
                    "Calamagrosis psoudo phragmatis" = c(90,    102,    63),
                    "Carex diluta" = c(109, 64, 82),
                    "Carex divisa" = c(96,  99, 59),
                    "Carex pseudocyperus" = c(130,  71, 54),
                    "Carex stenophylla" = c(97, 98, 59),
                    "Catabrosa aquatica" = c(100,   5,  150),
                    "Centaurea iberica" = c(124,    85, 46),
                    "Cirsium hygrophilum" = c(158,  42, 55),
                    "Cladium mariscus" = c(159, 96, 0),
                    "cod2" = c(54,  82, 119),
                    "Cynodon dactylon" = c(121, 54, 80),
                    "Eleocharis palustri" = c(124,  100,    31),
                    "Epilobium parviflorum" = c(67, 80, 107),
                    "Eromopoa persica" = c(83,  15, 157),
                    "Funaria cf.microstoma" = c(8,  0,  247),
                    "Glaux maritime" = c(4, 196,    55),
                    "Hordeum brevisubulatum" = c(76,    70, 109),
                    "Hordeum glaucum" = c(40,   79, 136),
                    "Inula britannica" = c(95,  108,    51),
                    "Juncus articulatus" = c(107,   79, 69),
                    "Blysmus compressus" = c(81,    127,    47),
                    "Juncusinflexus"= c(149,    106,    0),
                    "Medicago polymorpha" = c(60,   86, 109),
                    "Mentha spicata" = c(150,   23, 82),
                    "Ononis spinosa" = c(66,    112,    77),
                    "Phragmites australis" = c(234, 0,  21),
                    "Plantago amplexicaulis" = c(108,   83, 64),
                    "Poa trivialis" = c(90, 28, 138),
                    "Polygonum paronychioides" = c(20,  12, 223),
                    "Potentila reptans" = c(106,    41, 108),
                    "Potentilla anserina" = c(105,  58, 91),
                    "Ranunculus grandiflorus" = c(129,  25, 101),
                    "Schoenus nigricans" = c(143,   91, 21),
                    "Setaria viridis" = c(10,   7,  238),
                    "Sonchus oleraceus" = c(178,    0,  77),
                    "Taraxacum officinale" = c(117, 28, 110),
                    "Trifolium repens" = c(94,  4,  157),
                    "Triglochin martima" = c(63,    96, 95),
                    "Veronica anagallis-aquatica" = c(55,   37, 163)
)
AddToTernary(points, data_points, pch=21, cex=1.2,
             bg=vapply(data_points,
                       function (x) rgb(x[1], x[2], x[3], 128, maxColorValue=255),
                       character(1))
)

AddToTernary(text, data_points, names(data_points), cex=0.8, font=1)

```

# Answer
> 1 votes
Hmmm I'm not familiar with AddtoTernary function, I can see that you've placed a function within the "AddtoTernary" function. The traditional ggplot2 way of doing this would simply be controlled by the theme

```
ggtern(data=x,aes(x1,x2,x3)) +
  geom_point(fill="red",shape=21,size=3) +
  theme_tern_bw()

```

This will give an arrow for each axis, but additional code is required for different colours (I'd have to look up my past work).

I think what you might want is a **colour key**? In this case ggtern will have a "key" function because ggplot2 has one.Here, your colouring is based on species not on a colour gradient as to how far to the corner of a graph you are.

Alternatively you can easily put names besides each point in ggplot2, but again this is traditional style code.

Also you can add different shapes to the graph, which might be useful instead of colours.

Arrows and names together peronsally I've not seen on a ggplot2 graphic and I'd have not idea how to do it.

FYI I don't think we need all the data just a few points fine, although it shows why you prefer arrows - personally I would use colour coding. In "traditional" ggplot2 coding the graphing is built up in layers, so if additional functionality is required it is simply appended to the exisiting command. Thus it is quite easy to make modifications.

---
Tags: r, ggplot2
---


# Question
Title: What is "aligned sequences" and "consensus sequence" in the context of sequence logo? How to compute these?
In bioinformatics, a sequence logo is a graphical representation of the sequence conservation of nucleotides (in a strand of DNA/RNA) or amino acids (in protein sequences).

so, What is "aligned sequences" and "consensus sequence"?

how can I get these sequences based on following data?

Figure 2.5 from the book "Machine Learning - A Probabilistic Perspective"

more info about this figure at here

# Answer
Panel (a) shows a sequence gap-free alignment. Each row corresponds to a contiguous 15 base pair sequence of DNA (e.g. row 1 could be a human sequence, row 2 could be the equivalent mouse sequence etc.). Each column corresponds to a particular position/residue (e.g. column 1 is the 1st position of the alignment etc.).

Panel (b) shows the sequence logo. As noted in the comments, this does not correspond to the alignment in (a). Instead, here is the corresponding logo:

From this you can intuitively understand that at the 1st sequence position you have only As and Ts and their frequency corresponds to the y-axis. For the 3rd position, you have only As, so only A is represented here in the logo etc.

The consensus sequence is the most frequent letter at each position in the alignment, i.e. `[AT]TAGC` etc. where `[AT]` means A or T since they are jointly high in frequency (5 each).
> 5 votes

---
Tags: sequence-alignment, ngs, sequence-analysis
---


# Question
Title: Why do ten rows (Figure_1) correspond to 2 bits (Figure_2) in a sequence logo?
Following this question, I'm confused with the computation of sequence logo

Following data comes from the book "Machine Learning - A Probabilistic Perspective (Figure\_1)"

here is the corresponding sequence logo (Figure\_2).

ten rows represent sequences of DNA (e.g. row 1 could be a human sequence, row 2 could be the equivalent mouse sequence etc.)

Each column corresponds to a particular position.

Why do ten rows (Figure\_1) correspond to 2 bits (Figure\_2)?

# Answer
> 6 votes
Bits are not frequencies. If a position only contains an A (position 3 for example) then you would need 2 questions (bits) to derive that information without a priori knowledge.

1. Is it a G or C? if no \> then it is a A or T
2. Is it T? If no then it is an A

Position 1 can be derived by 1 question only:

1. is it a G or C? If no then it is an A or T

In this case the number of bits can be calculated using the frequenci)es.

**For 4 bases, A only (100%)**:

4 * -0.25 * log2(0.25) - (-1 * log2(1)) = 2 (bits)

**For 8 letters, A only (100%)**:

8 * -1/8 * log(1/8)- (-1 * log2(1)) = 3 (bits)

**Mixed base possibility (A/T 50% each) with 4 bases**:

4 * -0.25 * log2(0.25) - (-0.5 * log2(0.5)+-0.5 * log2(0.5)) = 1 (bit)

As mentioned in the comments, you can follow the wiki link for an example.

Or, an example with a little more explanation from the biology SE site.

sequence logo calcuation

Hopes this helps

---
Tags: sequence-alignment, ngs, sequence-analysis
---


# Question
Title: Problem: "Pair-end" reads scRNA seq data (Drop-seq)
In case of Drop-seq, we have paired end data.

Read 1: Cell code + UMI (unique molecule identifier)

Read 2: The transcript information

But I have a problem/doubt with the sample I am working on.

The sample I am using is the following (Check the "Reads" tab).

As you know the Drop-seq is "paired-end", we are expected to see two reads for a spot. Although this sample say paired-end, it has only one read per spot.

For example I can share a link of a different scRNA-seq data where you can properly see two reads for a spot

Example sample, you need to check the "Reads" tab.

Where I am going wrong?

I asked one of the main authors of the paper. The following is the reply I got :

> I recommend that you download the aligned BAM files that are hosted in the same GEO record. Read 1 is already processed into the cell and UMI barcodes and held as custom tags (XC and XM) in the BAM files. The cells are already barcode-corrected, so if you use those files, your cell barcodes will line up with mine; if you start from FASTQs, they will not. For most aligners, you can just use the BAM file as input to realign. (It has all reads, even those that did not align.)

But I could not find any "XM" or "XC" keywords in the bam file :(

To understand his reply you have to be familiar the processing steps of the Drop-seq.

Looks like they have submitted some kind of processed data. I could not figure out how much the data is processed. I am trying to use the data starting at different processing steps but I could not figure out how much the data is processed.

# Answer
> 3 votes
I received the following message from the author of another Drop-seq paper:

> After much conversation with the folks at the SRA, it appears the SRA compression algorithm automatically deletes certain subsets of bam tags and does so in cryptic way (in the sense of not alerting the uploader). Obviously this is going to be an issue for many single cell datasets encoded this way, and I hope the SRA will eventually change their system to accommodate.

If you look at the "Download" tab of the same page, there now exists a link to download the originally uploaded BAM file.

Alternately, you can use the SRA AWS tutorial to download the files on that platform.

---
Tags: scrnaseq, read-mapping
---


# Question
Title: How to compute the Shannon entropy for a strand of DNA?
I'm confused by the computation of sequence logo. Wikipedia gives a process about this without a concrete example.

Let's just consider DNA, so there are 4 bases (nucleic acids).

The following data comes from the book "Machine Learning - A Probabilistic Perspective (Figure 1)"

The Shannon entropy of position $i$ is:

${\displaystyle H\_{i}=-\sum \_{b=a}^{t}f\_{b,i}\times \log \_{2}f\_{b,i}}$

Where ${\displaystyle f\_{b,i}}$ is the relative frequency of base

This post is computing position 3, where it seems that the relative frequency of `a` is 100%, and the relative frequency of 3 other bases are 0%.

Now, how do you calculate the Shannon entropy of position 3? It seems to be 0, which is apparently incorrect, where am I wrong?

# Answer
Why do you think the entropy of 0 is incorrect? It intuitively makes sense, as there is no uncertainty about the base at position 3, and thus there is no entropy.

However, what is plotted in a sequence logo isn't the entropy, but rather a measure of the "decrease in uncertainty" as the sequence is aligned. This is calculated by taking the entropy at this position if we randomly aligned sequences ($H\_g(L)$), and subtracting from it the entropy of the alignment ($H\_g(s)$): $$ R\_{sequence}=H\_g(L) - H\_s(L) $$

The entropy at position 3 based on a random alignment is calculated by assuming there are 4 equally likely events (one for each base), and thus: $$ \begin{align\*} H\_g(L) & = -((1/4 \times -2) + (1/4 \times -2) + (1/4 \times -2) + (1/4 \times -2)) \\ H\_g(L) & = 2 \end{align\*} $$ Notably, $H\_g(L)$ will always be 2 when dealing with nucleotide sequences.

So in your example, we have: $$ \begin{align\*} R\_{sequence}&=H\_g(L) - H\_s(L) \\ R\_{sequence}&=2 - 0 \\ R\_{sequence}&=2 \text{ bits of entropy} \end{align\*} $$

Finally, to work out the height of each base at each position in the logo, we multiply the frequency of that base by the overall information gain $R\_{sequence}$. In this case the frequency of the base A at position 3 is of course 1: $$ \begin {align\*} H(b,l) &=f(b, l) \times R\_{sequence} \\ H(A, 3) &= f(a, 3) \times 2 \\ &= 1 \times 2 \\ &= 2 \end{align\*} $$

Source: Sequence logos: a new way to display consensus sequences. (Thomas D. Schneider, R. Michael Stephens, 1990)

Note: I've excluded the correction factor $e(n)$ from the overall information-gain calculation for the sake of simplicity. Refer to the paper above for an explanation of this.
> 5 votes

---
Tags: sequence-alignment, ngs, sequence-analysis
---


title: Subset a multisample VCF file

text:
# Question
Title: Subset a multisample VCF file
I need some help with GATK. I have a multisample vcf file containing 18125 samples. After performing PCA I have around 303 samples as outliers and I wanted to remove these samples from my multisample vcf file.

I tried doing it using the GATK `SelectVariant` tool. It does some stuff and generates a VCF file which is a bit smaller in size than the original file, but in the end the total count of samples remains the same. As I have no prior knowledge of GATK (working with it for the first time) I am just wondering if you guys can help me out in fixing this issue.

The command that I used was

```
gatk SelectVariants -V filtered.vcf -xl-sn remove.txt -O test.vcf

```

`remove.txt` contains the sample names in each row (1 sample per row).

After running this command it generated a `test.vcf` file, for which I counted samples using:

```
bcftools query -l test.vcf | wc -l

```

However, I got exactly same number of samples as in the original VCF file `filtered.vcf`.

# Answer
Based on the way you've called the command (using `-xl-sn`), I assume you are using GATK 4 and not GATK 3.

If this is the case, the documentation for the `-xl-sn` flag is available here: https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org\_broadinstitute\_hellbender\_tools\_walkers\_variantutils\_SelectVariants.php#--exclude-sample-name.

By default, the `-xl-sn SAMPLE` means "exclude the sample called SAMPLE from my VCF". In your case, you're providing an entire list of samples, in a file, which is not the default behaviour. In order to tell GATK that you're providing a file and not a sample name, **you need to ensure the file you provide has the extension `.args`**.

From the documentation:

> File names must use the extension ".args"
> 3 votes

---
Tags: vcf, quality-control, gatk
---


# Question
Title: How to retrieve the subcellular location info from uniprot?
I have a list of proteins with their gene names/symbols. I would like to identify which of them are extracellular proteins or plasma membrane proteins.

I can get this information in UniprotKB by searching the proteins one by one. Because I have hundreds of proteins, I wonder whether I can do a systematic search of the UniprotKB in RStudio. I tried the UniProt.ws select method. But I couldn't even get the example code in the package work. I downloaded the uniprot\_sprot.dat file. But I don't know how to read this file correctly in RStudio. I used tab-delimited format. It turned out there is only one variable in the file.

Can anybody tell me what I should do to solve the problem? I am pretty new to any computer languages.

# Answer
> 2 votes
Since you are already on UniProtKB, the simplest way will be if you add this information to the output table and download it as a that you can import and filter in R.

To get additional annotations, click the button labelled "Columns" on top of the UniPortKB website.
This allows you to add a wide variety of annotations such as subcellular location to the output. Tick the checkmark for subcellular location and save on the top right.
Then you can download the proteome-wide annotation for example as a tab-separated list using the "Download" dialog next to the "Columns" selection. This list can be imported into R.

# Answer
> 0 votes
After my post of the question, I went to try the UniProt.ws package in the RStudio again. I realized that there is a mistake in the example code in the package. The up object is from a human database, but the keys the example provides are the ENTRZID from mouse proteins. When I tried human ENTRZID numbers, it worked. The only problem for the select method in the UniProt.ws is that it can't retrieve the protein information for a large batch of proteins. The maximum number that worked for me is about 80 proteins per time.

---
Tags: r, sequence-annotation, uniprot
---


# Question
Title: what does a sequence look like before alignment?
I am confused with Sequence alignment, which is

> a way of arranging the sequences of DNA, RNA, or protein to identify regions of similarity that may be a consequence of functional, structural, or evolutionary relationships between the sequences.

to understand this easily, I would like to use this image (image\_1)

the image above (image\_1) is the primary part of nucleic acid structure (primary, secondary, tertiary, and quaternary), shown below (image\_2)

is the sequence gcgc... an aligned sequence? if yes, what does the sequence look like before alignment?

# Answer
> 4 votes
Sequence alignment typically refers to the primary sequence. That is, we ignore the complex 3D structure of the molecule and its spatial arrangement in the cell. Instead, we use a symbolic representation composed of symbols like A, C, G, T, and U. Alignment is a conceptual tool that, as you pointed out, is used to find regions of similarity between sequences. Therefore, it only makes sense in the context of two or more similar sequences.

By way of example, consider a simple case of two sequences X and Y. This is what they would look like before alignment. The similarities between the two sequences are unclear.

```
X: ATTGACATAGGCACTTCGACTA
Y: GACATCTTCGGCTAA

```

This is what they would look like after alignment. The similarities (and differences) are now much easier to see.

```
X: ATTGACATAGGCACTTCCACTA-
      |||||     ||||  |||
Y: ---GACAT-----CTTCGGCTAA

```

This arrangement identifies regions of similarity that *may* be related through evolution. This doesn't always mean perfect agreement though. As this example shows, sometimes mismatches or gaps are allowed so that long stretches of identical matches are lined up correctly. Sometimes, for very closely related sequences, the optimal alignment doesn't require any gaps at all. This seems to be the case with some of your recent questions about sequence logos.

**Note**: Sequence alignment is a fundamental concept in bioinformatics, and there is a tremendous amount of research and development on the topic. It's hard to imagine doing any amount of serious DNA sequence analysis responsibly without a basic understanding of sequence alignment. I would recommend searching Google and/or Google Scholar for intro-level reviews of the topic.

---
Tags: sequence-alignment, ngs
---


title: Multiple line chart

text:
# Question
Title: Multiple line chart
I am trying to create a multiple coloured line graph based on six lines. Essentially I have data that identifies different units of vegetation by the variable” NDVI (vegetation index)” I want to create a plot with differently coloured lines for each units.

The following is a screen shot of the data structure and my code. It is plotting all of the data. But the time data on X-axis is disrupted and not ordered. What am I missing? And also how do I get legend for my chart?

```
library(ggplot2)
library(tidyverse)

df2 <- data.frame(Time<-c("28Jan","13Feb","16Mar","1Apr","17Apr","3May","19May","4Jun","20Jun","6Jul","22Jul","7Aug","23Aug","8Sep","24Sep","10Oct","26Oct","11Nov","13Dec","29Dec"),
                  A<-c(0.209,0.225,0.237,0.312,0.377,0.468,0.499,0.492,0.45,0.43,0.405,0.405,0.378,0.358,0.337,0.297,0.179,0.305,0.241,0.183),
                  B<-c(-0.472,-0.589,-0.422,-0.479,-0.371,-0.112,-0.03,-0.134,-0.132,-0.12,-0.1,-0.227,-0.137,-0.025,-0.02,-0.011,-0.024,-0.377,-0.494,-0.49),
                  C<-c(-0.278,-0.385,-0.31,-0.3,-0.279,-0.033,0.089,0.154,0.231,0.308,0.341,0.4,0.397, 0.383,0.37,0.309,0.248, 0.271,-0.212,-0.3),
                  D<-c(-0.544,-0.604,-0.505,-0.541,-0.475,-0.173,-0.054,-0.026,0.098,0.199,0.257,0.338,0.37,0.372,0.36,0.274,0.217,-0.007,-0.552,-0.528),
                  E<-c( 0.213,0.239,0.277,0.304,0.384,0.47,0.537,0.547,0.5,0.489,0.452,0.469,0.439,0.419,0.408,0.372,0.049,0.335,0.304,0.282),
                  F<-c(0.234,0.248,0.3,0.324,0.391,0.444,0.512,0.503,0.448,0.434,0.407,0.421,0.4,0.385,0.376,0.335,0.087,0.323,0.295,0.275))
head(df2)

ggplot(df2, aes(x = Time, group = 1)) +
  geom_line(aes(y = A), colour="blue") +
  geom_line(aes(y = B), colour = "red") +
  geom_line(aes(y = C), colour = "green") +
  geom_line(aes(y = D), colour = "black") +
  geom_line(aes(y = E), colour = "yellow") +
  geom_line(aes(y = F), colour = "brown") +
  ylab(label="NDVI") + xlab("Time")

```

# Answer
Couple of things regarding your code:

* It would be good to have your data frame in long format for ggplot2
* For having the x-axis ordered chronologically, you have to specify it manually, otherwise it is ordered alphabetically.

Then, here is my solution:

```
library(ggplot2)
library(reshape2)

# I have substituted your <- by = and changed the name of the groups
df2 <- data.frame(Time = c("28Jan","13Feb","16Mar","1Apr","17Apr","3May","19May","4Jun","20Jun","6Jul","22Jul","7Aug","23Aug","8Sep","24Sep","10Oct","26Oct","11Nov","13Dec","29Dec"),
                  g_A = c(0.209,0.225,0.237,0.312,0.377,0.468,0.499,0.492,0.45,0.43,0.405,0.405,0.378,0.358,0.337,0.297,0.179,0.305,0.241,0.183),
                  g_B = c(-0.472,-0.589,-0.422,-0.479,-0.371,-0.112,-0.03,-0.134,-0.132,-0.12,-0.1,-0.227,-0.137,-0.025,-0.02,-0.011,-0.024,-0.377,-0.494,-0.49),
                  g_C = c(-0.278,-0.385,-0.31,-0.3,-0.279,-0.033,0.089,0.154,0.231,0.308,0.341,0.4,0.397, 0.383,0.37,0.309,0.248, 0.271,-0.212,-0.3),
                  g_D = c(-0.544,-0.604,-0.505,-0.541,-0.475,-0.173,-0.054,-0.026,0.098,0.199,0.257,0.338,0.37,0.372,0.36,0.274,0.217,-0.007,-0.552,-0.528),
                  g_E = c( 0.213,0.239,0.277,0.304,0.384,0.47,0.537,0.547,0.5,0.489,0.452,0.469,0.439,0.419,0.408,0.372,0.049,0.335,0.304,0.282),
                  g_F = c(0.234,0.248,0.3,0.324,0.391,0.444,0.512,0.503,0.448,0.434,0.407,0.421,0.4,0.385,0.376,0.335,0.087,0.323,0.295,0.275))

# Convert to the proper format
df2.melted <- melt(df2)

# Order of the factor in x-axis
df2.melted$Time <- factor(df2.melted$Time, levels = c("28Jan","13Feb","16Mar","1Apr","17Apr","3May","19May","4Jun","20Jun","6Jul","22Jul","7Aug","23Aug","8Sep","24Sep","10Oct","26Oct","11Nov","13Dec","29Dec"))

# Note group = variable to connect the points of the different groups
ggplot(df2.melted, aes(x = Time, y = value, colour = variable, group = variable)) +
  geom_line() +
  ylab("NDVI") + xlab("Time")

```
> 3 votes

---
Tags: r, ggplot2
---


# Question
Title: UCSC liftover cannot produce bed format output
I was trying to use UCSC liftover tool to convert human coordinates to mouse. While most of the files worked, one file didn't:

```
chr16   86430087    86430726    Human_chr16:86430087-86430726___element_1___positive____neural_tube[12/12]___hindbrain_(rhombencephalon)[12/12]___limb[3/12]___cranial_nerve[8/12]
chr16   85620095    85621736    Human_chr16:85620095-85621736___element_2___negative_
chr16   80423343    80424652    Human_chr16:80423343-80424652___element_3___negative_
chr16   80372593    80373755    Human_chr16:80372593-80373755___element_4___positive____neural_tube[6/10]___hindbrain_(rhombencephalon)[10/10]___midbrain_(mesencephalon)[10/10]

```

Does anyone know why this bed file can only be converted to chr:start-end format but not bed4 format?

# Answer
> 1 votes
I figured out the issue by myself. It turned out that UCSC liftover doesn't like the character ":"

I removed ":" and it worked

---
Tags: bed, ucsc, liftover
---


title: Intersecting multiple data frames

text:
# Question
Title: Intersecting multiple data frames
I have to intersect multiple data frame to find common elements. For visualization purpose i can use upset Plot but it doesn;t return a object containing the intersection object .I tried jvenn but at max it can take six sets to intersect .

Im able to do the interesction but i would like to know if i can do it more efficiently ,which apparently not sure how to do this is my code as of now where i read each file as data frame, covert them to list and pass them to do intersection .

```
NewList <- split(GO0000122, f = seq(nrow(GO0000122)))
NewList1 <- split(GO0006351,f=seq(nrow(GO0006351)))
NewList2 <- split(GO0006355,f=seq(nrow(GO0006355)))
NewList3 <- split(GO0006357,f=seq(nrow(GO0006357)))
NewList4 <- split(GO0006366,f=seq(nrow(GO0006366)))
NewList5 <- split(GO0030154,f=seq(nrow(GO0030154)))
NewList6 <- split(GO0045892,f=seq(nrow(GO0045892)))
NewList8 <- split(GO0045893,f=seq(nrow(GO0045893)))
NewList7 <- split(GO0045944,f=seq(nrow(GO0045944)))

a <- intersect_all(NewList,NewList1,NewList2,NewList3,NewList4,NewList5,NewList6,NewList7,NewList8)

b <- do.call(rbind.data.frame, a)

```

Can i read the files and pass them into list and do the same that would help me from doing lot of things manually

I can read the files from the directory where i have my files but not sure how to pass them to make it as list.

```
filenames <- gsub("\\.csv$","", list.files(pattern="\\.csv$"))
for(i in filenames){
  assign(i, read.csv(paste(i, ".csv", sep="")))
}

```

The object filenames is `character` here .It just returns me list of files.

So how to make my question bit more clear , i would like to read the list of file from the directory and then pass these to perform intersection. Any suggestion or help would be really appreciated .

# Answer
Change dataframe to load all the files into a list:

`wd <- "C:/all_your_dataframes_in_one_folder"`

`setwd(wd)`

Then you get the list of files you want from your folder

`files <- list.files()`

Put into a list with a loop

```
dflist <- list()
for(i in 1:length(files)){
 dflist[[i]] <- read.csv(files[i])
}

```

Then use tidyverse solution

`library(tidyverse)`

```
map(dflist, ~.$GeneID) %>% #creates a list just of the column of interest
    reduce(intersect)  #applies the intersect function cumulatively to the list

```
> 1 votes

---
Tags: r
---


# Question
Title: How to convert a .fastq chip-seq file into a wig/bigwig track and how to extract chip-seq peaks?
I have a reference .fasta file and a raw .fastq file with chip-seq data. I am trying to create a bigwig track from the and .fastq and .fasta ref file of the raw signal. Then I would like to do some peak analysis for the chip-seq track.

Can someone help me with googling the write things, I feel really lost while trying to google it and not sure where to start from. I was thinking to do `.fastq -> .bam -> wig -> bigwig` and then somehow to do peak detection analysis.

# Answer
Before you proceed, it is usual to do quality control on the .fastq data. The FASTQC program is a good place to start.

Once you have checked the quality, trimming is the typical next step to remove low-quality data. There are options such as Trimomatic, trimGalore for this.

Once you are sure you have clean data, your first step is alignment - align the sequences in the .fastq file to the reference sequence in the .fasta file. bowtie2 is one option for this process. This gets you to the .bam stage.

These steps will take a little time, some research, and may be somewhat frustrating, but they are the foundation for what you want to do. One option for somewhat automating these steps is to use the Galaxy platform (usegalaxy.org), otherwise these are mostly command-line tools.

The bam to wig and wig to bigwig steps are relatively simple conversions, and Biostars has answers for these.
> 1 votes

---
Tags: fastq, chip-seq
---


title: Identifying Indels from Chromatograms

text:
# Question
Title: Identifying Indels from Chromatograms
I have around 100 chromatograms (`.ab1` files) from Sanger sequencing a genome at loci believed to have an indel.

I'm new to interpreting this kind of data in general, but I've read a bit on the general idea—mostly in guides like this. What I'm not able to figure out from the resources I've checked is **how to identify insertions and deletions in a chromatogram**.

I'm most concerned with heterozygous indels, and these seem to be simpler. (And since my samples are (supposed to be) larger indels that are more rare and more likely to be het, I think.) A heterozygous indel, according to my understanding, would like something like consecutive heterozygous SNPs—overlapping peaks, maybe not with the same height, but with the same horizontal placement—until the end of the sample.

In other words, something like the top track in the open window in this image.

(Source: CodonCode.)

(I'm trying to do this manually, and not with packages like this for several reasons, including: I want to see how this is actually done, I want to do it as thoroughly as possible—I will probably check my results with an automated package—and I was asked to do it like this.)

*How can evidence of (het) indels like this be distinguished from noise in the chromatogram?* My guess is it includes things like

* the overlapping peaks last until the end of the sample (and don't start at the beginning, but somewhere in the middle); and
* the peaks are still at the same horizontal positions; i.e., one is not shifted right or left relative to the other?

Is this correct?

---

With my own chromatograms, I have a lot of doubt about distinguishing noise from indel. Here's an example, visualized with 4Peaks. (Full chromatogram here.)

On the left, the chromatogram appears high quality. There are sharp, regularly spaced, single peaks. On the right, there appear to be two sets of peaks lined up with each other. Does this constitute evidence for an indel? Does the homopolymer region—I've noticed in many of these samples, the putative indel seems to start in a homopolymer region—cast doubt on this?

---

I've also noticed that in some situations like this, there appear to be not two but *three* peaks.

Is one of these just noise and can be disregarded? Or does this cast doubt on the validity of the data—which is supposed to be coming from the diploid human genome?

# Answer
If you want to call heterozygous SNPs and INDELs, you may try polyphred. I believe there are a couple of other similar tools, but I have forgotten their names.

> How can evidence of (het) indels like this be distinguished from noise in the chromatogram?

Your two points are correct. Basically, at some position, you start to see two streaks of overlapping peaks. Once identified, the two streaks can be aligned to each other if you shift one of them. This is a strong constraint, but it is hard to tell by eye because it is not obvious how to assign peaks to streaks without running an automatic tool. Doing this is even more difficult given noisy traces like your last figure. I think you'd better run polyphred or similar.

> the putative indel seems to start in a homopolymer region—cast doubt on this?

The INDEL mutation rate in long homopolymer runs is orders of magnitude higher than elsewhere in the genome. This has two implications: 1) you expect to see more homopolymer indels; 2) PCR polymerase is more likely to introduce errors in homopolymers––a higher fraction of homopolymer indels tend to be wrong.
> 4 votes

# Answer
to understand how to interpret sanger sequencing chromatograms I'd like to give a brief introduction to how this data are generated (sorry, if you already know this)

The target region is amplified by PCR. The following sequencing reaction is quite similar to a PCR, despite two important differences:

1. Only one strand is used as a template, by adding just 1 primer instead of two.
2. In addition two normal dNTPs, we have modified ddNTPS. Those are fluorescence labeled and the 3'-OH is missing, so there is no chance to add more dNTPS.

The result of this reaction are single stranded DNA fragments of different length. The last base is fluorescence labeled. So we have fragments that looks like this:

```
A
AT
ATC
ATCG
ATCGA

```

Of course there a many copies of these fragments.

Now these fragments are separated by capillary electrophoresis. Small fragments are detected first. We can differ the several bases by its fluorescence color.

Let's imagine a DNA strand like this:

```
ATCGA

```

What happens if we had on the second strand a SNP?

```
ATCGA
ATAGA

```

Now the sequence would detect two signal of equal intensity on base 3 - one for C and one for A. This have two effects in the chromatogram.

1. Two peaks of equal heights with different colors
2. The relative peak height is just half the height, compared to those sequences that have just one base at this position.

Let's see what happen in case of a heterozygous deletion. An alignment would look like this

```
ATCGA
AT-GA

```

But our capillary electrophoresis doesn't know about alignment. It only knows about fragment length. So this is what it will detect

```
ATCGA
ATGA

```

With the beginning of the deletion it will detect to different signals. Just like in case of the SNP, the relativ signal intensity is now just about the half.

With this knowledge let's go back to your initial questions:

> How can evidence of (het) indels like this be distinguished from noise in the chromatogram? My guess is it includes things like
>
> * the overlapping peaks last until the end of the sample (and don't start at the beginning, but somewhere in the middle); and

They don't start at the beginning, if your indel is not right after the primer, correct. In case of an insertion or deletion you have strands of different size. So if your target region is quite small you may have luck, that you can see the end of the overlap at the end of the sequence, because the shorter strand reaches its end.

> * the peaks are still at the same horizontal positions; i.e., one is not shifted right or left relative to the other?

Yes. And one of the best things to distinguish between noise and true heterozygous variant is a relativ drop off of the signal intensities by approximately 50% compared to homozygous samples.

Concerning your second chromatogram: This is not an indel. As the other already said, this is a mistake made by the polymerase. Invitro polymerases doesn't like homopolymer stretches. They begin to add extra bases or "forget" one. As we do PCR this error accumulates. In my experience up to 9 time the same base is not a big problem, 10 only if you have luck, from 11 upwards your sequence are terrible. If you are interested in the sequence part behind this stretch, you have to sequence from the other direction.

fin swimmer
> 3 votes

# Answer
Sanger sequencing produces a bit of messy signal at the start and a lot at the end of the sequence, but should have a region of clean signal in the middle. If the signal is getting mixed up in the middle, it suggests that there are sequences that are out of phase with each other. That *could* be a heterozygous INDEL, but it could also indicate that multiple template sequences have been amplified by the same primer.

Sanger sequencing is frequently carried out in both directions (forward and reverse complement), and if similar things are happening in both directions at the same place, it makes it more likely that there is a heterozygous INDEL and not something else funny going on. It's a very good idea to have homozygous templates for all expected variants, if possible. If the secondary peaks in the signal match what is observed in the homozygous variants, that can also be used as evidence. In the absence of control sequences, it is also possible to eyeball the peaks to find self-similarity.

Where triple peaks are seen, that could point to a region where there are two INDEL variants in close proximity, or multiple regions getting amplified at the same time, or even a single INDEL variant with more than two variant states (e.g. a homopolymer with 8 As that is sometimes 7 As, and sometimes 9 As).
> 1 votes

# Answer
Honestly I would not try to call heterozygous indels with Sanger sequencing. It's just not the right technology.

If you were trying to confirm a known indel that didn't flank a homopolymer stretch, and you had sequence from the other side that got messy at exactly the same point, maybe you could safely conclude there was something going on there, and possibly you could eyeball the exact right sequence of mixed peaks, but I don't see that what's what you have in your examples at all. It just looks like noise.
> 1 votes

# Answer
There are specialized tools to deconvolute (separate) heterozygous InDels in a chromatogram trace using a reference genome. For instance, Indigo is a web application for this, the command-line application is called Tracy.
> 1 votes

# Answer
One important point that @gringer briefly touched on is that there is often off-target amplification when performing PCR on whole-genome DNA.

So, the **noise in the chromatograms very well could be from another part of the genome and not related to indels at all**.

The only way to ensure that this is not a problem is to clone your PCR products into a vector, transform into something like E. coli, grow up a single colony, cut the product with restriction enzymes, size-select, then sequence the correct-sized insert. Or you could just re-sequence your PCR product with Illumina or PacBio to accurately sequence the complete population of molecules -- indels and all. Check out:

Hebert, Paul DN, et al. "A Sequel to Sanger: amplicon sequencing that scales." BMC genomics 19.1 (2018): 219.

The safe thing to do is to not read too much into the noise of Sanger chromatograms unless it is a simple and very clean SNP.
> 0 votes

---
Tags: sequence-alignment, ngs, genome, variant-calling
---


# Question
Title: Is there a convention or standard of the highest bits in a sequence logo?
In bioinformatics, a sequence logo is a graphical representation of the sequence conservation of nucleotides (in a strand of DNA/RNA) or amino acids (in protein sequences).

wiki says

> The total height of the letters depicts the information content of the position, in bits.

and then gives this example

> A sequence logo showing the most conserved bases around the initiation codon from all human mRNAs (Kozak consensus sequence). Note that the initiation codon is not drawn to scale, otherwise the letters AUG would each have a height of 2 bits.

it seems that there is a 0.1 bit letter.

is this in the context of 2bits? is 2bits the highest bits in this sequence logo? Is there a convention or standard of the highest bits in a sequence logo?

# Answer
> 1 votes
This is explained in the original paper \[1\] describing the logos (emphasis mine):

> The importance of a particular position in a binding site is more clearly and consistently given by the information required to describe the pattern there \[7, 8\]. To choose one symbol or state from two equally likely possibilities requires one "bit" of information. For example, to communicate the result of a coin-ip to someone requires 1 bit of information because only one yes-no question needs to be answered: "Is it heads?". **If a position in a binding site always contains one base (e.g. G), then we need exactly two bits of information since two yes-no questions need to be answered: "Is it A or G?" (i.e. is it a purine?) and "Is it A or C?". (If the answers to both questions are "no", it must be a T.)** Furthermore, if a position contains two bases (e.g. sometimes A and sometimes G), only one question suffices since a two out of four choice is equivalent to a one out of two choice. Therefore, only one bit is needed to describe a position in a binding site that contains only purines, but two bits are needed to describe a position that always contains adenine.

So yes, 2bits is the highest value you will see on the y-axis of a ***DNA*** sequence logo (proteins are a different story since there are more possibilities). In the one you show, the scale is from 0 to 0.2, but it's the same principle.

However, you should think of logos as a useful visual aid, and not a way of getting rigorous mathematical information about a sequence. So don't put too much weight on the bit values. If you need that sort of data, calculate it from the alignment directly.

---

\[1\] Schneider TD, Stephens RM., Sequence logos: a new way to display consensus sequences., *Nucleic Acids Res.* 1990 Oct 25;18(**20**):6097-100.

---
Tags: ngs, sequence-analysis
---


# Question
Title: Subsetting of object existing of two samples
I have a Seurat object existing of an aggregate of two samples namely; RD1 and RD2. I am trying to make a subset of each sample. But it's not really working. I manage to make a subset of the cells but I don't manage to make a subset of the whole object. Isn't there a way to do this, based on the subset from cell names a subset to make from the whole object.

```
seuratObj$orignal_1 <- seuratObj@meta.data$orig.ident
seuratObj$orignal_2 <- seuratObj@meta.data$orig.ident
seuratObj$orignal_1 <- gsub("RD1", "TRUE",seuratObj$orignal_1 )
seuratObj$orignal_1  <- gsub("RD2", "False",seuratObj$orignal_1  )
seuratObj$orignal_2  <- gsub("RD1k", "TRUE",seuratObj$orignal_2  )
seuratObj$orignal_2  <- gsub("RD2", "False",seuratObj$orignal_2  )

subset_DR1 <- data.frame(orignal_1 = seuratObj$orignal_1)
subset_DR1 <- subset(subset_DR1, orignal_1 == TRUE,select = c("orignal_1"))

subset_DR2 <- data.frame(orignal_2 = seuratObj$orignal_2)
subset_DR2 <- subset(subset_DR2 , orignal_2 == TRUE,select = c("orignal_2"))

SubsetData(seuratObj, cells = rownames(subset_DR1), subset.name = seuratObj$orignal_1)

```

I just started to work with seurat packages. I think I am missing something.

# Answer
I think your issue is simply not assigning the newly-subset object in to a new object.

Assuming `rownames(subset_DR1)` is indeed a list of cell names

```
seuratObj_subset_dr1 <- SubsetData(seuratObj, cells=rownames(subset_DR1))

```

will work.
> 2 votes

---
Tags: r, seurat, subset
---


title: Clustering using cellranger

text:
# Question
Title: Clustering using cellranger
Is it possible to run cell ranger clustering using only a gene x cell matrix? Cellranger reanalyze performs clustering but requires a hd5f file. I don't want to run the whole pipeline (including alignment). I was trying to manually generate a hd5f file to input to cellranger reanalyze. Is this possible?

# Answer
The direct answer is yes. If you have never worked with hdf5 you can start here.

But, cellranger outputs several hdf5 files. You can extract the count matrix from the "/outs/filtered\_gene\_bc\_matrices\_h5.h5" file.

They provide a lot of information on how to use the hdf5 format so their their documentation is a good place to start.

In addition: you can also load the data into Seurat, monocle or scanpy and follow their tutorials to analyze the data.
> 2 votes

---
Tags: clustering, cellranger
---


# Question
Title: How to reproduce this specific sequence logo from 10 000 human mRNAs?
this is A sequence logo showing the most conserved bases around the initiation codon from 10000 human mRNAs.

I would like to reproduce this sequence logo. How can I get the sequences? which database should I use?

this search is apparently wrong, since I don't know the right way.

https://www.ncbi.nlm.nih.gov/nuccore/?term=10000+human+mRNAs

# Answer
If you look at the Appendix from https://www.ncbi.nlm.nih.gov/pmc/articles/PMC306349/, there is a table of sequences that were likely used to construct the sequence logo, at least partially, from the 5' end to the stereotyped ATG subsequence.

It's not in a nice format amenable to use directly with weblogo (like FASTA), but with some programming/elbow grease it might be rendered into that form, which you can then use.

Generally, you might investigate the supplementary materials of these and related papers, to see if you might something directly usable.
> 1 votes

---
Tags: sequence-alignment, ngs, sequence-analysis
---


# Question
Title: What is the name of this type of figure?
This figure comes from this wiki page.

I googled "gene location figure" and "gene location plot", none gets similar results. I also tried search by image, none of the results is close.

so, what is the name of this type of figure?

# Answer
> 6 votes
This kind of visualization showing the banding patterns of chromosomes are called "chromosome ideograms".

You can use for example IdeoViz in R to generate them.

# Answer
> 3 votes
Such a figure is called ideogram, or chromosome ideogram.

---
Tags: ngs, visualization, genomics, figure-explanation, genetics
---


# Question
Title: Access workdir defined on command line from within Snakefile
Snakemake provides access to a `workflow` object within a Snakefile. This allows one to, for example, have dynamic programmatic access to the directory containing the Snakefile (via the `workflow.basedir` attribute). Is there a similar way to access the working directory that is specified on the command line using `--directory`?

Through a bit of exploration I was able to figure out that the `workflow` object has a `workdir()` method, but it appears this is used internally by Snakemake for registering working directories, not for access.

**UPDATE**: It looks like the working directory is stored internally as the `_workdir` attribute of the `workflow` object. However, this attribute is unset (default value `None`) when I'm creating static/global objects in the Snakefile. Is it impossible to access the working directory prior to building the workflow DAG?

# Answer
> 2 votes
You could query Python for the working directory within the Snakefile, since evidently Snakemake changes the actual working directory of the process to the one specified by `--directory`. For example, using a stub Snakefile with just these two lines:

```
import os
print("Current working directory: " + os.getcwd())

```

And on the command line:

```
$ snakemake --directory working-directory
Creating specified working directory working-directory.
Current working directory: /home/jesse/dev/snakemake-wd/working-directory
Building DAG of jobs...
Nothing to be done.
Complete log: /home/jesse/dev/snakemake-wd/working-directory/.snakemake/log/2019-08-02T142009.349338.snakemake.log

```

---
Tags: python, snakemake, workflow-management
---


# Question
Title: Huge variability in blastn speed between same short sequence job
I have been running a couple of blast searches and I am very suprised and somewhat frustrated by the huge variability in job completion time. The specifics of my blast search are as follows:

* Running blast using webbrowser (chrome) via ncbi
* 19 nucleotide query
* database: WGS of one particular sequencing project/bacterial strain (so database is VERY small, in total approx 4.5 MB divided over 60 contigs)
* parameters are mostly untouched and are automatically optimized for short queries by blast

In one instance, I started the job more than an hour ago and it is still not done. Knowing from experience that this previously hasn't taken that much time, I started a new job in a different browser tab with exactly the same sequence, database and parameters, and now it finished within 20 seconds. What on earth is happening here?

# Answer
First of all, if you're using an online service like NCBI's BLAST, then the time to completion will depend on the load the service is experiencing. Your job is ut in a queue and will be started when it's its turn. If the service is under heavy load, then its turn might come later than you expect. Then, if the service is under heavy load, even after launching your job, it might take longer to finish because it won't have enough resources to run quickly.

Since you say that the same job finished quickly on a new tab, it sounds like either the first job was last in the queue or there was just a bug that caused it not be launched or something like that. We can't really help with things that resolve themselves like that.
> 1 votes

---
Tags: blast
---


title: Error Viewing file on HiGlass

text:
# Question
Title: Error Viewing file on HiGlass
I am implementing hi-glass python on Jupyter, and upon running the test data code

```
from higlass.client import View, Track
from higlass.tilesets import cooler

ts1 = cooler('../data/Dixon2012-J1-NcoI-R1-filtered.100kb.multires.cool')
tr1 = Track('heatmap', tileset=ts1)
view1 = View([tr1])
display, server, viewconf = higlass.display([view1])

display

```

I get the error message stating:

```
Error Retrieving tileset info: http://localhost:60501/spi/v1

```

Please advise on the correcting. The code above has been copied from the HiGlass webpage.

# Answer
I'm one of the authors of HiGlass. Sorry you're getting that error. I just tried it out in a fresh environment and it worked for me. Just one thing is missing in the docs which is that you need to do `import higlass`.

Here's what I did to set up my environment:

```
conda create --name test-bi python=3.6
source activate test-bi

```

The rest is from the docs:

```
pip install jupyter higlass-python
jupyter nbextension install --py --sys-prefix --symlink higlass
jupyter nbextension enable --py --sys-prefix higlass
jupyter notebook

```

It looks like either you or somebody else added a bug report here: https://github.com/higlass/higlass-python/issues/29

I'll respond and continue the conversation there.

Also, there's a repository full of examples from the HiGlass SciPy presentation here: https://github.com/higlass/scipy19
> 1 votes

---
Tags: visualization, genomics
---


title: Biopython Failure to break out of loop

text:
# Question
Title: Biopython Failure to break out of loop
After qblast, I wanted to print only details from the first 3 BLAST output. Yet, the code below printed the details of every BLAST output. Why can't I break out of the loop to achieve my desired output?

```
for alignment in blast_records.alignments:
    countResults = 0
    for hsp in alignment.hsps:
        if countResults == maxResults:
            countResults = 0
        break
        else:
            print(alignment.title)
            print(alignment.length/ len(seq)*100)
            print(hsp.expect)
            print(hsp.identities)
            countResults += 1

```

# Answer
> 2 votes
Assuming that your actual code is this:

```
for alignment in blast_records.alignments:
    countResults = 0
    for hsp in alignment.hsps:
        if countResults == maxResults:
            countResults = 0
            break
        else:
            print(alignment.title)
            print(alignment.length/ len(seq)*100)
            print(hsp.expect)
            print(hsp.identities)
            countResults += 1

```

Here you have an outer loop, which is iterating over all the possible database matches, and an inner loop, iterating over the alignments within that database entry. You are currently limiting the inner loop, but not the outer loop, which you probably want to do. Try this:

```
# Take the first 3 database hits
for alignment in blast_records.alignments[:3]:
    for hsp in alignment.hsps:
            print(alignment.title)
            print(alignment.length/ len(seq)*100)
            print(hsp.expect)
            print(hsp.identities)

```

# Answer
> 2 votes
There was an indentation bug (its a Python thing). Personally I'd use,

```
for alignment in blast_records.alignments:
    countResults = 0
    maxResults = 3
    for hsp in alignment.hsps:
        if not countResults == maxResults:
            print(alignment.title)
            print(alignment.length/ len(seq)*100)
            print(hsp.expect)
            print(hsp.identities)
            countResults += 1
        else:
            countResults = 0 # I'd dedinitely remove this in your original code
            break

```

The loop gives a bit more control if there is any subsequent filtering, species, percentage identity etc.. There are other ways to write it.

Here's an alternative

```
for alignment in blast_records.alignments:
    maxResults = 3
    for countResults, hsp in enumerate(alignment.hsps):
        while countResults <= maxResults:
            print(\
            alignment.title, "\n",\
            alignment.length/ len(seq)*100, "\n",\
            hsp.expect, "\n", hsp.identities, "\n"\
            )

```

BTW the enumerate function keeps track of the counting.

Alternative print ...

```
           out = """ alignment.title,
           alignment.length/ len(seq)*100,
           hsp.expect,
           hsp.identities """
           print (out)

```

Let me know if there are bugs

---
Tags: python
---