text:
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?

Note that the frequency is always in the interval [0,1] instead of [0,100]. H is 0 for columns with only one symbol.

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.

Thanks for your answer. What does "g" mean/represent/come from in $H_g(L)$

I'm not sure. My guess is g=genome and s=sequence(ing) but that's pure conjecture since the paper doesn't spell it out

Thanks for your reply. And does L mean/represent/come from location since the paper says "at various positions (L)".

Yes, that would be my assumption. This is important to include in the equation because it's only calculating the information at a single position. The authors propose another equation for an entire sequence (which is just summing these per-position scores)

Thanks a lot. And last question, does R come from/represent/mean acids?

I don't know. But also I think you're putting too much importance on the authors' choice of variables. You could write any of these equations with different letters and they would be just as correct.


text:
Store and search on variant and annotation data
I need to store a large number of variants with their annotations. I need to import dbSNP, as well as Clinvar, and ExAC etc. Also I have sample variant data from patients. I need to be able to search on genomic regions as well as annotations. And it needs to be scalable.
What would be a preferred data storage technology for such a case? What is the industry using? It needs to be accessible by APIs/REST.
We are using ElasticSearch, but that is proving to be expensive in terms of servers needed.
Something similar to this has been discussed in
The state, limitations and comparisons of large variant stores, but that post is a little old, I am sure there has been progress since then.

The question you linked is about storing genotypes. Your question is about keeping site information, which is different and much easier.
As to solutions, you can parse the information you want and store them in a standard SQL database. You can use UCSC binning for region queries. Modern SQL databases all provide R-tree indices, which can be adapted for region query and is probably better these days.

A standard SQL database will be slow for such large amounts of data. I will be filtering on annotations as well as regions.

@user2689782 Have you tried? SQLite for example easily handles tens of millions of variants at my hand. UCSC/EnsEMBL keep the entire dbSNP in MySQL as well. The most efficient way to deal with annotations highly depends on the type of annotations and queries. It is more about schema design.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
How to demultiplex a mix of single-indexed and dual-indexed samples
The problem
If I have a sample sheet that contains both single-indexed and dual-indexed samples, I can split it up into two sample sheets and then run bcl2fastq on each one. However, when doing this, large Undetermined fastq files are generated. E.g., when processing the single-indexed samples, all the dual-indexed samples go to Undetermined. And when processing the dual-indexed samples, the single-indexed samples go to Undetermined.
Additionally, because they are being processed separately, if any single index A is part of a dual index A+B, then when processing the single indexes, an A+B may be mistaken for an A, so it would seem that they need to processed simultaneously to avoid this mis-assignment.
The question
Given a sample sheet and directories of BCL files, how can such a set of sequencing data be demultiplexed correctly, either using bcl2fastq or the Picard tools?
To put it another way, I want to demultiplex a single sequencing run that contains both single-indexed and dual-indexed samples. It can be assumed that the indexes are sufficiently distinct such that any sample's index configuration is different from any other. But assuming that the different index configurations are not segregated to particular lanes of the sequencer, the question is how to demultiplex the files correctly such that both the single-indexed and dual-indexed samples are recognized.
Attempts at a solution
Using bcl2fastq directly
If a sample sheet of a mix of indexes is given to bcl2fastq (v2.17) directly, it produces the error
ERROR: bcl2fastq::common::Exception: Success (0): .../bcl2fastq2/src/cxx/lib/layout/BarcodeCollisionDetector.cpp(127):
Throw in function void bcl2fastq::layout::BarcodeCollisionDetector::validateNewBarcodeSizesAgainstExisting(const std::vector<long unsigned int>&) const
Dynamic exception type: boost::exception_detail::clone_impl<bcl2fastq::layout::BarcodeCollisionError>
std::exception::what: Barcodes have an unequal number of components.

Barcodes have an unequal number of components.

Using picard
It seems that it should be possible using the Picard tools, but I have not found a way to set up the inputs to ExtractIlluminaBarcodes and IlluminaBasecallsToFastq that processes this configuration correctly.
The syntax for supporting multiple index configurations is not entirely clear. But when using various combinations of N and '*' on multiplex_params.tsv and barcodes.txt required by the Picard tools, a large Undetermined file is still produced and actual sample fastq files are tiny, indicating it is not processing them correctly.
Using multiple calls to bcl2fastq
As indicated in the discussion above, this suffers from the problem of having all the dual indexed samples going to the "Undetermined" file when processing the single-indexed samples, or vice versa, and creates two output directories of reports and stats, which must be merged.
Padding empty indexes with N
By padding, I mean converting ACGTACGT to ACGTACGT+NNNNNNNN so that the single-index samples in the sample sheet are "dual" as well. This strategy seems that it would work except there is a bug in bcl2fastq that it treats "N" literally instead of as a wildcard. See the release notes for details.

Welcome to the site. Sorry but I don't understand which files do you have. Could you clarify what are your files? Or do you have a mix of single-indexed and dual-indexed and you can't differentiate between them?

@llrs Hi, thanks for the question. In this case, I'm trying to demultiplex a single sequencing run that contains both single-indexed and dual-indexed samples. You can assume that the indexes are sufficiently distinct such that any sample's index configuration is different from any other. But assuming that the different index configurations are not segregated to particular lanes of the sequencer, the question is how to demultiplex the files correctly such that both the single-indexed and dual-indexed samples are recognized.

Ok, now I understood the question, but I'm sorry I can't help with it. Perhaps you could [edit] your last comment in the question for other users.

I'm late to the party but simply run cellranger mkfastq twice, with different arguments for the mask and with --filter-dual-index for the double-indexed samples. In a second step you then have to disentangle the dual A+B reads from single-index A reads based on the readID's that are found in both of the demultiplexings. This is cumbersome and time+diskspace consuming, but it can be done. I would use some Unix tools for that but there are probably smarter ways.
EDIT: I mistakenly said bcl2fastq, this must be cellranger mkfastq. (thanks AmadeusDrZaius)

Thanks for your suggestion. I only see that "--filter-dual-index" option available for cellranger mkfastq, rather than bcl2fastq itself. Is that where you saw the option?

ah yes, that is true, sorry. I guess you would have to do an additional filtering step on the FASTQ file, the second word on the ID lines is 'index1+index2' for properly demultiplexed reads.

Have you just tried giving bcl2fastq one sample sheet with a mix of single and dual indices?  I don't think what you are trying to do is a problem.

if any single index A is part of a dual index A+B,

Well, it might be a problem if you did that.  That was poor planning.

That was the first thing tried but unfortunately bcl2fastq will give the error "Barcodes have an unequal number of components" when you try to do that (using v2.17). I'll add a note to my question.

Are you sure that the single index samples don't have a site for the second index to sit down and read?  My lab does some samples with single indicies and some with double, but they always prep with libraries with kits meant for dual indexing, so the second index is always there, even if I ignore it.

That's a good idea that might be a good workaround in some cases, but I think the general answer is that they are not always physically dual-indexed. The use of "N" as described in the bcl2fastq documentation would work if not for the fact that a bug in bcl2fastq causes it to treat "N" as a literal base "N" instead of as a wildcard. Given that, I'm guessing that the Picard tools are more likely to have a full solution for this particular situation.

I don't think it's going to be trivial for software to figure out the difference between "second index doesn't match what I expect" and "second index isn't really here".  You might have to have bcl2fastq export the indices as fastqs, so you can use the quality score to parse yourself.

That's an interesting point. I expected that it would just make a "best guess", but I can see why it wouldn't. It sounds like that would primarily only be an issue in the case of "A" and "A+B" that I described above. If the case of "A" and "A+B" is avoided, it should be able to tell the difference. Thanks for your thoughts.

I would be surprised if the Experiment Manager software would allow a run configuration that was un-processable.  I take it the instrument was run without a sample sheet made by Experiment Manager?

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
How to reverse complement the DNA sequences for given inverse/reverse coordinates?
I have the series of coordinates in id.txt file, whose coordinates sequences are in genome.fasta file.
The coordinates of id.txt file is shown below,
Contig3:15-7
Contig2:5-10
Contig1:12-3

The genome.fasta file is shown below,
>Contig1
AAGGCCATCAAGGACGTGGATGAGGTCGTCAAG
>Contig2
ACGAAGCGCTCGCCAAGGCCGAAGAAGAAGGCC
>Contig3
GCTGCGGCGCTGATCCTGGCGGCCCGCGCCGAG

I have used the following codes to extract the sequences from genome.fasta based on the coordinates in id.txt by using the following command xargs samtools faidx genome.fasta <id.txt > result.fasta
It has extracted only Contig2:5-10sequence. Because, the coordinates of Contig2:5-10is in proper order. But Contig3:15-7 and Contig1:12-3 coordinates are in inverse/reverse order, So samtools could not fetch those sequences.
I need to extract those sequences from Contig1:12-3 and Contig3:15-7 coordinates and also I want to reverse complement them. I have  larger fasta file and plenty of coordinates sequences to be extracted, So I need to automate this process. So, please help me to automate this process.
The expected output is
>Contig3
GATCAGCGC
>Contig2
AGCGCT
>Contig1
CTTGATGGCC

Thank you in advance.

Welcome to the site. Do you have any preference about the tool to be used? What else have you tried to set the coordinates in the right order? Usually this kind of problems are solved by working on every tiny problem one after the other.

I have used linux command the following command to keep the coordinates in right order,                                                                                "awk '$2 > $3 {printf( "%s\t%s\t%s\n", $1, $3, $2); next} id.csv".

So now you have the right order and the positions. What is stopping you from extracting the sequences?

By using the above mentioned command I formatted the coordinates in right order and made it compatible for samtools. But still, one more thing I have to do is reverse complement of those right ordered sequences. I have jumbled coordinates, so the command changes the order wherever it needed.

I need to extract the sequence based on the coordinates specified above (even though it is mixed up with inverse ordered and right ordered). But, Now I need to reverse complement the right ordered (right ordered by using above command) sequences. So that I could go for my further analysis. Moreover, I need to automate this process. The reverse complement implementation to be taken place, wherever the inverse coordinates are there.

What do you know biopython? It was designed to make reverse-complement sequences and translate if needed

The exonerate tool comes with some very handy utilities for sequence manipulation. Here, we are interested in fastasubseq and fasterevcomp. First, install exonerate as described here. If you're using Debian or Ubuntu or one of their derivatives, it might be enough to run
sudo apt install exonerate

Then, follow these steps:

Create one fasta file for each sequence in genome.fasta
fastaexplode genome.fasta

A little preprocessing
$ awk -F[:-] '{
                if($2>$3){
                    start=$3-1;
                    len=$2-start" -"
                }
                else{
                    start=$2-1;
                    len=$3-start
                }
                print $1,start,len}' id.txt
Contig3 7 8 -
Contig2 5 5
Contig1 3 9 -

So, that will print the target contig's name, the start coordinate of the target sequence, the length of the target sequence and a - if it is on the reverse strand. The fastasubseq program I will use below starts counting from 0 and not from 1, so we also need to subtract 1 from the coordinates to get them right.
Now, feed that into fastafetch and fastarevcomp
tmpFile=$(mktemp);
    > subseqs.fa
    awk -F[:-] '{
                if($2>$3){
                    start=$3-1;
                    len=$2-start" -"
                }
                else{
                    start=$2-1;
                    len=$3-start
                }
                print $1,start,len}' id.txt |
    while read cont start len rev; do
        fastasubseq "$cont".fa $start $len > $tmpFile;
        if [[ -n $rev ]]; then
            fastarevcomp "$tmpFile" >> subseqs.fa;
        else
            cat "$tmpFile" >> subseqs.fa;
    fi
    done

This will produce the following file when run on your data:
$ cat subseqs.fa
>Contig3:subseq(6,9) [revcomp]
GATCAGCGC
>Contig2:subseq(4,6)
AGCGCT
>Contig1:subseq(2,10) [revcomp]
CTTGATGGCC

Finally, remove the per-contig fasta files
rm Contig*fa

Thank you so much terdon for your time and help. Your script works awesome. However, I did a mistake in my question. Instead of mentioning Contig3:15-7, Contig1:12-3, I have mentioned Contig3:15-7 and Contig2:5-10.
I need to reverse complement the sequences, only if it in reverse order/inverse order. Which ever sequences are not in reverse order (Properly order sequences), only to be extracted as its not to be reverse complemented.

Please accept my sincere apologies for the same. Could you please do the needful.

@K.Dineshkumar I'm sorry, but I don't understand the problem. The script will only reverse complement if the start position is smaller than the end position. Isn't that what you want? In the example above, Contig3:15-7 and Contig1:12-3 have been reverse complimented, but Contig2:5-10 has not.

@K.Dineshkumar please edit your question and add this. Show us the expected output. It's very hard to read in the comments.

@K.Dineshkumar
Never mind, I found the problem. The fastasubseq tool uses 0-based coordinates (the first position of a sequence is 0), but you have 1-based (the first position of a sequence is 1). Try the updated answer.

@terdon Please see my edited question, wherein I have included my expected outcome.

@terdon, The revised/updated script printing one base less in the end of the sequence. For example,  Contig2:5-10 says, print the sequence from 5 to 10 But the script print only 5-9.

@K.Dineshkumar argh! Yes, of course, sorry. Fixed now, try the (last, I hope) update. And thanks for adding the output to your question, now I can check :)

@terdon, Now it works perfect. I really appreciate your help. Thank you.

@K.Dineshkumar you're welcome. If this answers solves your issue, please take a moment to accept it by clicking on the checkmark on the left. That will mark the question as answered and is the way that thanks are conveyed on the Stack Exchange sites. You don't need to do it now, you might want to wait in case anyone can give a better solution, but please remember to accept an answer since that marks the question as answered.

@K.Dineshkumar thanks for accepting! Feel free to accept one of the other answers though, I won't take it the wrong way :) They give more efficient solutions.

Since you tagged python, here is a solution that solves this in a single script using Biopython:
import sys
from Bio import SeqIO

d = {}
with open('id.txt') as f:
    for line in f:
        contig, coords = line.split(':')
        start, end = map(int, coords.split('-'))
        d[contig] = start, end

for record in SeqIO.parse('genome.fasta', 'fasta'):
    start, end = d[record.id]
    if start > end:
        record.seq = record.seq[end - 1: start].reverse_complement()
    else:
        record.seq = record.seq[start - 1: end]
    SeqIO.write(record, sys.stdout, 'fasta')

Output on your example:
>Contig1
CTTGATGGCC
>Contig2
AGCGCT
>Contig3
GATCAGCGC

Thank you @Chris_rands

I try to use it for a bacterial genome. But it shows error as follows,
Traceback (most recent call last):   File "/home/ga/Documents/commands/scripts/amphora/reverse_comple_seq_extrac.py", line 12, in <module>     start, end = d[record.id] KeyError: '1'

@K.Dineshkumar try using d[record.description] instead, i can't know what will work without seeing your data though

My genome of interest can be obtained from the following link https://www.ncbi.nlm.nih.gov/Traces/wgs/CAHO01?display=contigs

My coordinate is as follows,                                                                                CAHO01000012.1:92353-95295 CAHO01000017.1:27362-27949 CAHO01000014.1:11449-12012 CAHO01000018.1:51156-52928 CAHO01000030.1:58079-56889 CAHO01000001.1:60348-62039 CAHO01000026.1:45206-44508

When I tried with d[record.description] instead of d[record.id], I end up the error as follows,                                                                                Traceback (most recent call last):   File "/home/ga/Documents/commands/scripts/amphora/reverse_comple_seq_extrac.py", line 12, in <module>     start, end = d[record.description] KeyError: 'CAHO01000001.1 Xanthomonas citri pv. mangiferaeindicae LMG 941 WGS project CAHO00000000 data, contig XMINcontig00001, whole genome shotgun sequence'

If you're okay with using a Perl script, you could convert your headers into stranded BED6 and use the following script bed2faidxsta.pl with samtools-indexed FASTA to get back stranded sequence:
https://gist.github.com/alexpreynolds/fa9b0f90e181e3b4f640
This script looks at the sixth column of your BED file of chromosome, start, and stop positions, to decide whether to do a reverse complement of the queried sequence.
You could use awk to quickly build such a BED6 file:
$ awk -vFS='[:-]' '($0~/^>/){ if ($2 < $3) { strand = "+"; } else if ($2 > $3) { strand = "-"; } else { printf("Coordinates must differ\n"); exit 1; } printf("%s\t%d\t%d\t.\t0\t%c\n", $1, $2, $3, strand; }' in.fa > out.bed6

To set up indexed FASTA, here's a generic shell-based approach, which uses [Kent utilities][1] fetchChromSizes to get a list of chromosomes for the build mm9:
$ mkdir mm9
$ cd mm9
$ while read chr; do echo "$chr"; wget -qO- "http://hgdownload.cse.ucsc.edu/goldenpath/mm9/chromosomes/$chr.fa.gz" | gunzip -c - > $chr.fa; samtools faidx $chr.fa; done < <(fetchChromSizes mm9 | cut -f1)

This can be adjusted for mm10, hg19, hg38, and so on.
Once you have indexed FASTA and the above Perl script, you could do something as easy as the following:
$ ./bed2faidxsta.pl --fastaDir="/path/to/indexed/mm9" < in.bed > out.fa

If you're working with a reference genome, the nice thing about this approach is that awk, Perl, and indexed lookups are fast. You can reuse indexed FASTA for other projects that also require sequence lookups by genomic interval.

Thank you @Alex Reynolds


text:
Issue with visualising cladogram/phylogenetic tree with multiple sequence alignment data in R?
I would like to visualize tree with multiple sequence alignment.
My din.newick file is shown below,
(org1:0.49999997,org2:0.47298786,org3:28.37858179);

My din.fasta alignment file shown below,
>org1
ATGCTA
>org2
TAGCAT
>org3
AGTAGT

I have tried the following code,
library(tidyverse)
library(ggtree)
library(seqinr)
tree = read.tree("din.newick")
msaplot(p=ggtree(tree), fasta = "din.fasta")
msaplot(p=ggtree(tree) + geom_tiplab(align=TRUE), fasta = "din.fasta")

The above mentioned code is generating tree with the alignment matrix, but the tree root tip labels are not visible. It is either overlapped with matrix or it is not at all visible.  Kindly help me to fix this issue.

You were very close, you just need to supply the offset argument to msaplot().
So using your code block above with one extra argument in msaplot()
library(tidyverse)
library(ggtree)
library(seqinr)
tree = read.tree("din.newick")
msaplot(p=ggtree(tree), fasta = "din.fasta")
msaplot(p=ggtree(tree) + geom_tiplab(align=TRUE), fasta = "din.fasta", offset = 5)

yields

Thank you for supplying a minimal and reproducible example!

Thank you @Kohl Kinning.

I suspect the solution is to realign the two plots within "ggarrange", available here
The problem is stacking ggtree with geom_tiplab. The ggarrange library is described on STHDA here
I have not used this library however, but it should be easy because you are working singly within R. I tend to use ggplot2 within Python via rpy2, which is much more difficult and restricted.


text:
How to obtain the extended unaligned bases of the read in reference based genome assembly?
I did illumina read mapping based on a particular gene sequence (reference sequence) using the following command,
bowtie2-build -f ref_gene.fasta ref_index && bowtie2 -p 4 -N 1 -t -x ref_index -1 genome_R1.fastq -2 genome_R2.fastq -S tal.sam && samtools view -bS tal.sam | samtools sort - -o tal.bam && samtools mpileup -uf ref_gene.fasta tal.bam | bcftools call -c | vcfutils.pl vcf2fq > cns.fastq && seqtk seq -aQ64 -q20 -n N cns.fastq > contig.fasta

It maps the reads according to the reference sequence and yields, exactly same length of aligned/assembled sequences. But, I need obtain the unaligned overlapping read bases at the both end (reverse and forward end)  Please see the image for better understanding.

Is there any way to get the same?

You're not going to get vastly different replies here than Biostars or github.

@DevonRyan sniper!


text:
Pattern mining from a genomic sequence
I need to find the following pattern from a genomic sequence
TTCGB-N15-TTCGB
B represents C, G, or T
N15 represents A, T, G or C and these four base combination should have the length of 15 bases.

To do the same, I tried grep command, but I could not specify B and N15 in the grep command line. Therefore, please help me to find the same.
For example, I have a nucleotide sequences as follows
>ATGCTTATTGCCCATTTCGTGCATGCATATGCGCATTCGCGATCGATTAGGGATAT

Expected outcome,

TTCGTGCATGCATATGCGCATTCGC

Thanks in advance.

Could you please update your question to clarify that N15 represents A, T, G or C for exactly 15 bases?

echo ">ATGCTTATTGCCCATTTCGTGCATGCATATGCGCATTCGCGATCGATTAGGGATAT" | grep -oP  'TTCG[CGT][ATGC]{1,15}TTCG[CGT]'

[CGT] looks for one occurrence of one of the three.
[ATGC]{1,15} looks for "up to 15 letter long" combinations of the four bases.

Moreover, -P makes your regex call Perl-like and is required to make this regex work. Different regex flavors behave (slightly) differently.

Dear @haci,                                                                                I have followed your command line for the test data (I had given in the example section of my question), which gives proper result. But when I tried the same with my genome file with the following command, echo $(<org15.fasta) | grep -oP  'TTCG[CGT][ATGC]{1,15}TTCG[CGT]' and I am getting output as follows,   TTCGCCGCGGCGACGTTTTCGC TTCGCATTCGG TTCGCCACCTTGTTCGG TTCGGCTTTCGC TTCGTGGCCTTCGG TTCGGCTTCTTCGT TTCGTCAAGGACCCGATTTCGG TTCGCATCTTCGC

The thing is that, the script mining TTCGB.....TTCGB, but fails to fulfil the length N15 criteria. kindly help me to fix this. Thank you.

What is your "N15 criteria" though? [ATGC]{1,15} looks for up to 15 repetitions of the four bases as you specified in your question. And as far as I can see in the output you get, this is what you are getting, 2 TTCGB stretches are connected by up to 15 bases.

N15 usually means exactly 15 of any nucleotides. So use [ACGTN]{15} . Also, -i flag for grep may need to be added if some of the sequence is lower case, for example repeat masked.

@haci N15 means any of the 4 bases, but N15 length should be 15. exactly, The output should be like this, TTCGB-ATGCTAGCATGATAG-TTCGB. The output sequence should have TTCGB at both ends, at the same time in between TTCGB 15 bases also should be there. That 15 bases can be the combination of any of the four bases, but it should have 15 bases. Sorry, if my question is not clear.

OK, then just change [ATGC]{1,15} to [ATGC]{15} as Timur Shtatland suggested. The reason for [ATGC]{1,15} was your phrase in the question: "N represents A, T, G or C up to 15 bases".

@haci I did it echo $(<org15.fasta) | grep -oP 'TTCG[CGT][ATGCN]{15}TTCG[CGT]'. It works, thank you so much @haci and @Timur Shtatland.


text:
Proof of Breakpoint Reversal Sorting Approximation Algorithm
I'm kind of new to bioinformatics and trying to self-study. I'm reading a bioinformatics book: An Introduction to Bioinformatics Algorithms and ran into some problems about understanding the proof of the Breakpoint Reversal Sort algorithm in chapter 5.4.
Specifically, the book states a theorem:

If a permutation $\pi$ contains a decreasing strip, then there is a reversal $\rho$ that decreases the number of breakpoints in $\pi$, that is, breakpoint ($\rho(\pi)$) < breakpoint($\pi$).

This theorem guarantees that the algorithm will terminate and lead to no breakpoint. The book later explains the proof for the theorem:

Let $k$ be the smallest element of decreasing strip, so element $k-1$ would be the ending of an increasing strip. Therefore element $k$ and $k-1$ corresponds to a breakpoint. Then reversing the segment between $k$ and $k-1$ will decrease the number of breakpoints.

I'm quite confused about this proof and don't really understand what the authors mean by reversing the segment between $k$ and $k-1$. Can someone please explain the proof for me?
Here's the lecture slide about reversal sorting that based on the book http://csbio.unc.edu/mcmillan/Media/Lecture07Spring2015.pdf.

Check out this lecture for an easier understanding.

Denote the strips in your permutation as $\langle S_1, S_2, \cdots, S_n \rangle$. Suppose $k-1$ is in segment $S_i$, and $k$ is in segment $S_j$, for some $i, j \in [n], i \neq j$.
Recall by definition of $k$ that $k-1$ must be the maximum and therefore rightmost element in an increasing strip, and $k$ the minimum and therefore also rightmost in a decreasing strip. Consider the two cases:

$i < j$. The increasing strip comes before the decreasing. $k-1$ is on the 'right' end of $S_i$, but $k$ is on the 'right' end of $S_j$ also. In order for $k-1$ and $k$ to be adjacent, reverse the segments $S_{i+1}, \cdots, S_j$.
$j < i$. The decreasing strip comes before the increasing. $k$ is on the 'right' end of $S_j$, but $k-1$ is on the 'right' end of $S_i$ also. Again, we reverse $S_{i+1}, \cdots, S_j$.

So, "reversing the segment between $k$ and $k-1$" means to reverse all the strips in-between the two strips containing $k$ and $k-1$, excluding the 'leftmost' strip, but including the 'rightmost' strip.


text:
What do the symbols mean in minimap2's gap cost equation?
Heng Li's paper on minimap2 is understandable to me right up until the moment that it gives this definition of gap cost:

2.2.1 Alignment with 2-piece affine gap cost
Minimap2 performs DP-based global alignment between adjacent anchors in a chain. It uses a 2-piece affine gap cost (Gotoh, 1990):
$$\gamma_a(l) = \min \{ q + |l| \cdot e, \tilde{q} + |l| \cdot \tilde{e} \}$$

Suddenly, out of nowhere, we've got five variables used without prior introduction: $l$, $q$, $e$, $\tilde{q}$, and $\tilde{e}$.
Heng doesn't state their meanings, and based on a quick skim and text search of the cited paper by Gotoh, that paper doesn't use these variable names anywhere either.
What does this equation mean?

Here is the translation:

$l$ = gap length
$q$ = initial gap penalty for shorter gaps
$e$ = short gap extension penalty per every additional gap
$\tilde{q}$ = initial gap penalty for longer gaps
$\tilde{e}$ = long gap extension penalty per every additional gap

To understand the first of the two expressions being minimized, see the Wikipedia section on affine gap penalties.
Notice the similar notation structure
$$A + B \cdot L$$
where

$A$ is the known gap penalty,
$B$ is the gap extension penalty, and
$L$ is the length of the gap

To understand the other variables with the tildes on top, see the text below Equation 2 (emphasis added to key explanation):

On the condition that $e > \tilde{e}$,
  it applies cost
  $q + |l| \cdot e$
  to gaps shorter than
  $[(\tilde{q} - q)/(e - \tilde{e})]$
  and applies
  $\tilde{q} + |l| \cdot \tilde{e}$
  to longer gaps. This scheme helps to recover longer insertions and deletions (INDELs).

In other words, we make the assertion that we should penalize longer gaps differently than shorter gaps.
As a bonus, here's a visual understanding of why we use the minimization function on those two expressions using Python and some made-up (and probably not biologically useful) penalty terms. Note the requirement that $e > \tilde{e}$.
# Import libraries
import numpy as np
import matplotlib.pyplot as plt

# Setup
l = np.array(range(0, 15, 1))
y_short = eval('10 + 5*l') # Q=10, E=5
y_long = eval('20 + 3*l')  # Q=20, E=3

# Plot
plt.plot(l, y_short)
plt.plot(l, y_long)
plt.legend(["Short Gap Penalty", "Long Gap Penalty"])

Using our example and the expression $[(\tilde{q} - q)/(e - \tilde{e})]$, it will apply the first expression in the minimization function (the one without the tildes) when the gap length is shorter than 5.
$$
[(\tilde{q} - q)/(e - \tilde{e})]
=
[20 - 10)/(5 - 3)]
=
10 / 2
= 5
$$
We can see from our plot that the two lines intersect at 5, which we've verified through some arithmetic.


text:
How can I correctly get 3D scatter plot for micro-array pheno data in R [updated]?
I intend to see data points distribution (each reows) within labeled groups (different batchs such as 1,2,3, and so on) in 3D scatter plot, because I want to see the distribution of the data points and want to see how similar each group of data points in 3D space. To do so, I used scatterplot3d package from CRAN to get 3D scatter plot, didn't get correct plot for my data.
reproducible data:
Here is the reproducible data that simulated from the original dataset.
persons_df <- data.frame(ID= LETTERS[1:20], ages=sample(1:50,20, replace = FALSE),
                    batch=c(rep(1,5), rep(2,6),rep(3,4),rep(4,5)),
                    platform=c(rep("CATH",12), rep("ANN",8)))

basically, I want to render 3D plot for the data points that belongs to the same batch with the same color (1 is red, 2 is blue, 3 is green, 4 is yellow).
this was my attempt to render PCA plot
require(scatterplot3d)
par(mar=c(4,4,4,4), cex=1.0, cex.main=0.8, cex.axis=0.8)
scatterplot3d(persons_df$batch)

but above code doesn't work because first I should group data points that belongs to different batches then render them in 3D space. This is not very intuitive to me how to do it in R? Any idea to render possible 3D scatter plot for my data? any thought?
desired 3D plot:
I want 3D plot something like this:

how can I make this happen? any idea to get this done?

PCA doesn't make sense here. PCA is a dimensional scaling technique. When we use PCA, we are trying to represent (usually) highly dimensional data in two or three dimesnions. For example, in the example you linked prcomp(iris[, -5]) is performing PCA on sepal length, sepal width, petal length, and petal width. So we have 4 dimensions here.
In your reproducible example, you only have two dimensions to look at, Age and Batch. Number of dimensions aside, I think you are missing the useful data on which it would be useful to apply PCA. This is just metadata.

@KohlKinning agreed, I am going to rephrase my question. I need 3D scatter plot to see data points (a.k.a, number of rows that belongs to different batch). Do you have any idea?


text:
Allele Count and Allele Frequency in VCF files
I'm working in bioinformatics, but my computational skills far outstrip my knowledge of biology or genomics. So forgive the noobish question.
According to the VCF specification, the INFO column for each coordinate [CHROM, POS] can contain values for
AC : allele count in genotypes, for each ALT allele, in the same order as listed
AF : allele frequency for each ALT allele in the same order as listed: use this when estimated from primary data, not called genotypes
AN : total number of alleles in called genotypes
For example, I'm working with the the VCF file you can download from Kaviar.
I think I understand how to calculate those numbers from raw sequences from a pileup. If for a given coordinate [CHROM, POS], $N_0$ is the number of hits for the reference allele, $\textrm{AC}_i \equiv N_i$ is the number of hits for the $i^{th}$ variant (maybe only one), and AN is the total number of alleles counted: $\textrm{AN} \equiv \sum_0^n N_i$, then the allele frequency for the $i^{th}$ variant is
$$\textrm{AF}_i = \frac{\textrm{AC}_i}{\textrm{AN}} = \frac{N_i}{\sum_0^n N_i}$$
So the phrase "total number of alleles in the dataset" would mean "total number of hits to coordinate [CHROM, POS] for any allele. And by "hit" I mean an aligned sequence fragment that includes that coordinate.
But I may not be taking into account the caveat for AF: "use this when estimated from primary data, not called genotypes."
Am I correctly understanding AF, AC, and AN in VCF files?

No, that is not quite right. specifically, this part:

And by "hit" I mean an aligned sequence fragment that includes that coordinate.

To get data on how many reads contributed to calling a specific allele at a given position, one would use the DP field, for example, which isn't included in the Kaviar VCF.
The AC, AN, and AF fields in a VCF file are meant to be used to show the allele's frequency in the context of all of the samples used when making that VCF file.
Investigation process
This is one line from the Kaviar vcf file:
1       10230   rs200279319     AC      AA,A    .       .       AF=0.0002654,0.0048525;AC=7,128;AN=26378;END=10231

If we look at the README attached with the Kaviar vcf file, it says,

"Kaviar is a compilation of SNVs, indels, and comlex variants observed in humans, designed to facilitate testing for the novelty and frequency of observed variants."

How many humans though? If we go to the Kaviar website it says:

Kaviar contains 162 million SNV sites (including 25M not in dbSNP) and incorporates data from 35 projects encompassing 77,781 individuals (13.2K whole genome, 64.6K exome).

If we look at the Kaviar VCF header we see that the fields are defined:
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele Count">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in data sources">

If we look at these fields again, AF=0.0002654,0.0048525;AC=7,128;AN=26378, the AN=26378 tells us that for this locus, the Kaviar dataset contains information from 26378/2=13189 humans. That's the 13.2K whole genomes the Kaviar website was talking about. This site I'm looking at must be a non-coding region, because there is no exome data (protein-coding regions).
For this same site, AC=7,128 means that out of 13189 humans (26378 haplotypes), we saw the variants 7 and 128 times. Out of 26378 haplotypes, this is a frequency of AF=0.0002654,0.0048525.

Thank you so much for that clear and detailed answer!


text:
Hypergeometic test for the overlap of two set of genomic intervals (e.g. from Chip-seq data)
Suppose I have done two Chip-seq experiments. Now after the bias correction step such as MACS2, I get two set of genomic intervals (or peaks, as what is usually called) that corresponds to each of the experiment. I want to test the significance of the overlap of these two set genomic intervals that corresponds to each experiment.
The common test employed seemed to me is the Hypergeometic test. There's even a package specifically aim for this task: ChiPpeakAnno.
Recall from the Hypergeometric distrbition -- the analogy is that we have a bag of marbles with two different colors, say one blue and another red. A sample from such distribution consists of marbles drawn from this bag, and a probability measure for this distribution is the number of blue marbles in the sample.
Could someone explain what is the correspondence between the hypergeometric test and the significance overlap between two different set of genomic intervals?

Interesting dicussion, its about sampling error in minority populations and whether this sampling error can throw a spurious "they are different result", when in fact its a sampling bias.

There is no such thing as a hypergeometric test, at least in statistical textbooks. It's a fisher test based on hypergeometric distribution.
If it is chip-seq for the same target, i.e biological replicates, significance of overlap is not quite meaningful. You get more information by for example checking the correlation of the coverage between your replicates using deeptools or anything else.
If you are comparing between two different targets, for example, chip for two different histone marks, you can use the fisher test or chi-square this way.

Partition the genome in X bp bins, these are all your marbles
Identity those that overlap with peaks from target 1 and those that don't
Do 2. for target2
construct a contingency table using 2 and 3.

So you should end up with something like this (purely made up):
                      no_overlap_target2 overlap_target2
  no_overlap_target1               9799             112
  overlap_target1                    88               1

And if there is overlap, the bottom diagonal with be higher.. less on the off-diagonal, and you can test it with a fisher.test.

So in your table, top left is the sum of the number of intervals in target1 and target2 that are not overlapping each other? Also the off-diagonal is the the number of target 1 that is not overlapping with target2 and vice versa?

sorry.. not very clear, this table is for example, we partition the genome into 10000 bins.. 9799 don't have target 1 or target 2, off diagonal is found in target2 but not in target1

Ah, I see. Thanks! Now I think I know how to do.

Yes, I have two targets. I think I'll just go ahead and use Granges package in R to do this task.

Do you mind if you could explain why this is corresponds to the hypergeometric distribution in this setting? In particular, why does "drawing the blue ball" from the urn corresponds to X-bp bins that are not in target1 nor target2?

https://en.wikipedia.org/wiki/Hypergeometric_distribution, look under working example, you are looking for the association between drawn and color

Just to make sure, I took a quick look at chippeakanno, and it uses a hypergeometric test for gene enrichment, not necessarily overlap of two chip-seq results.
As you mentioned, you can imagine in your set of genes, you have those that overlap with your peaks, and those that do not. Hypergeometric distribution calculates the probability of getting a certain number of genes that overlap with your peaks (which is the same as sampling some balls and calculate the p-value of some balls with one color). The hypergeometric test calculates a p-value to see whether this number of overlapping genes is more than expected.
I also read on Wikipedia that the P-value of hypergeometric test times 2 is exactly the P-value for Fisher's exact test. Fisher's exact test tests whether your overlap is more than (or less than) expectation. Hypergeometric tests just feel like a one-tailed version of Fisher's exact tests.


text:
How to highlight specific cells in Seurat 2.4
I used Seurat 2.4 on our scRNA dataset to obtain the following tSNE plot.I was able to successfully extract cell IDs from the different clusters, and generate gene expression profiles. The analysis, and the biology makes sense. I have been however stuck in trying to highlight specific cells we are interested in using the Cell IDs (barcodes). For e.g, in this plot I am hoping to highlight specific cells (~ 120 of them) within Cluster 5 (blue cluster on the top). It can be either in featureplot mode or in this plot itself by an overlay, it doesn't matter. All I have to show are the 120 cells within the cluster.
For eg. if cluster 5 has cells atag-1 atgc-2 atat-3, cacc-4 cat-5... i want to recreate this plot using atag-1 atgc-2, cacc-4
I tried using do.identify = TRUE but that is a visual selection and can select regions, but is not useful to select specific cells. I was thinking if there is a way to subset and replace cluster 5 to include just the cells (~120 odd out of the total 200) but was unable to crack it. I tried by using cells.use but then realized that doesn't work either.Any ideas will be greatly appreciated, to get me out of this rut.
Thank you

Try to do this with PAGODA2 from Jean Fan: https://github.com/hms-dbmi/pagoda2

You can "highlight" any cells you'd like by creating a new column (or modifying an existing column) in the meta data slot within the Seurat object. From there you can modify the contents of that vector to your liking and plot the results.
It's not clear if you've been successful in obtaining the barcodes for the subpopulation to begin with, because that's a separate issue all together, and you need to clarify what attribute you are trying to subset the cells on because there are a number of ways to accomplish that. If you've already obtained the barcodes then I'll describe how you would plot the results below.
Let's assume you'd like to keep the original cluster IDs for all other cells. We'll start off by duplicating the cluster ID column (this isn't necessary, but it's good practice to preserve the original IDs).
n <- seurat_obj@meta.data # for convenience
n$new_column <- n$res.1.0

Now you need to modify the ids for cluster 5 to break them up into 2 populations.
n[rownames(n) %in% my_barcodes, "new_column"] <- "5b"

You can plot the new column using DimPlot with the argument group.by = "new_column", which will make the subpopulation stand out from the rest of the cells in cluster 5.


text:
How to extract column from a matrix matching the another file (sample file)
I have two files master file.txt and sample.txt
master.txt
Name    GTEX.1117F.3226.SM.5N9CT    GTEX.111FC.3126.SM.5GZZ2    GTEX.111FC.3326.SM.5GZYV    GTEX.1128S.2726.SM.5H12C    GTEX.1128S.2826.SM.5N9DI    GTEX.117XS.3026.SM.5N9CA
ENSG00000223972 1   0   0   1   0   0
ENSG00000227232 298 168 197 106 221 184
ENSG00000278267 0   1   1   0   0   0
ENSG00000243485 0   2   0   1   0   1
ENSG00000237613 0   0   0   1   0   1
ENSG00000268020 1   0   0   1   0   1
ENSG00000240361 1   3   0   4   2   1
ENSG00000186092 3   1   1   1   1   2

Sample file.txt
GTEX.1117F
GTEX.111FC
GTEX.111XS

Desired output
Name    GTEX.1117F.3226.SM.5N9CT    GTEX.111FC.3126.SM.5GZZ2    GTEX.111FC.3326.SM.5GZYV    GTEX.117XS.3026.SM.5N9CA
ENSG00000223972 1   0   0   0
ENSG00000227232 298 168 197 184
ENSG00000278267 0   1   1   0
ENSG00000243485 0   2   0   1
ENSG00000237613 0   0   0   1
ENSG00000268020 1   0   0   1
ENSG00000240361 1   3   0   1
ENSG00000186092 3   1   1   2
ENSG00000238009 2   2   1   1
ENSG00000233750 0   3   3   0
ENSG00000268903 103 44  76  24

I tried using the R-code form this link
Print specific columns in a matrix on the basis of sample id's in the header
However, I am not able to figure out how to extract columns matching only half part of ID (like for GTEX.111FC I want both GTEX.111FC.3126.SM.5GZZ2 and   GTEX.111FC.3326.SM.5GZYV)
Thank you

Assuming your "master" files looks like this:
df <- data.frame(Name = LETTERS[1:5],
                 Col1 = sample(1:10,5),
                 Col2 = sample(1:10,5),
                 Col11 = sample(1:10,5),
                 Col3 = sample(1:10,5))

  Name Col1 Col2 Col11 Col3
1    A    5    4     3   10
2    B    1    6     5    6
3    C    3    1     9    4
4    D    9    9     8    9
5    E    8   10     1    1

For a sequence of pattern you want to extract, you can use base r function such as grep:
match_pat = c("Col1","Col2")

m <- unlist(sapply(match_pat,function(x) grep(x,colnames(df))))

Col11 Col12  Col2
    2     4     3

And then, extract your columns of interest by doing:
> df[,m]
  Col1 Col11 Col2
1    5     3    4
2    1     5    6
3    3     9    1
4    9     8    9
5    8     1   10

In your example, it will be something like:
match_pat <- c("GTEX.1117F","GTEX.111FC","GTEX.111XS")
m <- unlist(sapply(match_pat,function(x) grep(x,colnames(df))))

master[,m]
```

Thank you for your answer. However, my actual sample list has more than 100 id's and matrix has 2000 rows. I want match_pat to read a list of id's? Is it possible?

@Megha where does your list of ID's come from ? a file ?

@Law: I have a GTEX brain sample file (which have id's in the header like  GTEX.111FC.3126.SM.5GZZ2. I want to extract the age-related samples(20-60 yrs) from this file. The age-related data file is a list of files which says all id's starting with GTEX.111FC are 20-30 year old and so on. The issue is the id's in the list is the half part of ID in the main master file making it difficult for me to extract it

And yes age id file is a separated file with the list of ID's

Hi, I did it with bash. I got the full id's from the partial match.
xargs -i grep ^{} matched.Brain.txt.2nd.txt < 2ndcol.20_30 > new.file

Glad you succeed to solve your issue with Bash. To answer your question, you can define match_pat to be a vector of the column containing id files like `match_pat <- as.vector(list_id_file$NameOfTheColumWithID). Then, you can extract column form your GTEX brain sample file.


text:
How to statistically test the strength and weakness of a correlated pair?
Subset of df before groupby looks like
Diabetes_Group  DM-1    DM-2    corr
Group1  56  56  1
Group1  56  58  0.611
Group1  56  61  0.448
Group1  56  63  0.3587
Group1  58  56  0.611
Group1  58  58  1
Group1  58  61  0.60113
Group1  58  63  0.534858
Group1  61  56  0.448
Group1  61  58  0.60113
Group1  61  61  1
Group1  61  63  0.622206
Group1  63  56  0.3587
Group1  63  58  0.534858
Group1  63  61  0.622206
Group1  63  63  1
Group2  78  78  1
Group2  78  85  0.622703
Group2  78  94  0.622244
Group2  78  23  0.687826
Group2  85  78  0.622703
Group2  85  85  1
Group2  85  94  0.697313
Group2  85  23  0.612089
Group2  94  78  0.622244
Group2  94  85  0.697313
Group2  94  94  1
Group2  94  23  0.603834
Group2  23  78  0.687826
Group2  23  85  0.612089
Group2  23  94  0.603834
Group2  23  23  1
Group3  27  27  1
Group3  27  80  0.118955
Group3  27  147 0.32038
Group3  27  9   0.335264
Group3  80  27  0.118955
Group3  80  80  1
Group3  80  147 0.430287
Group3  80  9   0.406426
Group3  147 27  0.32038
Group3  147 80  0.430287
Group3  147 147 1
Group3  147 9   0.546452
Group3  9   27  0.335264
Group3  9   80  0.406426
Group3  9   147 0.546452
Group3  9   9   1

What statistical test can be used to do comparative analysis between the different ID's (DM1 or DM2) to check which ID pair has weaker/stronger correlation within each diabetes_group?

Is there any statistical test to categorise the data into three categories:
(a) pairs of DM1 AND DM2 which are highly correlated than others.
(b) pairs that are not highly correlated
(c) For the category 1 pairs (highly correlated), is there any particular ID that outperforms than other ID's among different groups?

This comprehensive analytical outline for the stats and all the stats questions could be captured inside 3 analyses. Implementing this will obviously take a bit of time, but there should not be any hidden caveats (transformations etc ...).
Question 3
The statistical test you want is K-means clustering based on the correlation score where the number of predefined groups is 3. The labels ... you will need to combine the group with the
Label          Value
Group1-56-56   1
...

This is an absolute solution to the data and really cool use of K-means.
Question 3c
An easy solution is:

perform a new K-means with e.g. 10 groups and find the group with highest correlation values.

A better solution is:

hierarchical clustering would identify specific groups and for this sub-question is a good analysis. The input is the correlation matrix, the output is a tree. HOWEVER, interpreting this could be difficult if you're not trained in reading trees.

Question 2
A 2-way ANOVA would be a really good solution to this data for Q2. The test would be the columns were DM1 DM2 and the rows Group1, Group2, Group3. Interpreting it would be fun because you'd also have to grapple with a concept called "interaction" - which could be important here. However, the ability to address a whole raft of questions within one analysis is what 2-way ANOVA would permit. There is a slight issue of unequal group sizes 1, 2 and 3, because it would be better if Group 1,2 and 3 etc... were equal in size. However, the raw data looks to be very amenable to ANOVA. You must not use the correlation index in this test - just the raw values.

K-means clustering is in scipy and scikit-learn
Hierarchical clustering is also in scipy
For ANOVA
import statsmodels.api as sm
from statsmodels.formula.api import ols

There will be easy approaches to 2-way ANOVA e.g. SPSS and I think Excel (needs a plug in). I suspect that Hierarchical clustering will also be in a nice point and click package like SPSS (I don't use them, so I don't precisely know.)

Is it possible for you to explain this in terms of code or any example code

Hi @Megha, you have asked an entire statistical project within a single question, so I answered as best I could.What I propose is to investigate the facets (K-means is a good start) and have a go and ask a specific question if you're stuck. I've updated the post with the location of the specific packages. You don't need to stick with Python, point and click packages are perfectly fine because this is end analysis. Beyond this there is no further process in a data pipeline


text:
Is the genome a programming language (i.e. LISP)? Can we analyze it with computer science theory?
Question moved here from the biology stack exchange.
It looks like the genome is a kind of LISP language.
Operon structure of genome:
- https://github.com/philschatz/microbiology-book/blob/master/resources/OSC_Microbio_11_07_Operon.jpg
LISP function structure:
- http://support.ircam.fr/docs/om/om6-manual/res/listprefix.png
Nested genes:
- https://player.slideplayer.com/31/9782148/data/images/img11.jpg
Nested LISP cons cells:
- https://www.tutorialspoint.com/lisp/images/treestructure.jpg
The basic similarity I am seeing is of general structure: (func var0 var1 ... varn) and nested cons cells.
The correspondence with genome structure is the operon portion, which has:

operator = function
structural genes = variables
structural genes can be nested like cons cells

Is this a global structure for the genome, then a specific prediction is we can do something simple like counting open and close markers in the nested genes and they will be equal.
So, why would such a correspondence matter, if true?  The basic idea is that if the LISP structure is consistent, then we could maybe look for other programming constructs built on the structure, such as recursion/looping, maybe templating or even object oriented programming. Or just identify simple properties like balanced opening closing markers. The mere fact it has a LISP structure would be fascinating if true.  And since LISP is essentially an abstract syntax tree (AST), which all programming languages have, perhaps this structure I see in the genome is actually an AST, which implies an even more fundamental correspondence between the genome and programming languages.
More generally, we have a wealth of knowledge reverse engineering human programming languages.  If the genome is also a programming language similar in some respects to human programming languages, such as being a LISP variant, then we could apply this same reverse engineering knowledge base to better understand the genome.
So, first of all, is my high level observation correct?  Does the genome, at least for bacteria, have a global LISP like structure?  If yes, then if the genome is a kind of programming language, i.e. LISP, can we analyze it like software to understand its functionality?  Has anyone done this before?  Simple searching on Google Scholar turns up lots of software for analyzing genomes, but nothing about analyzing the genome as software.  It would be interesting to know whether anyone has seriously examined the question before, and either made some useful discoveries, or found out there is only a very loose analogy with programming languages, and thus knowledge regarding programming languages is not of use regarding the genome.
Hopefully this question is clear enough now, if not please ask me for further clarification.  I would like to get some sort of decent answer.
UPDATE: Someone who has put thought into "DNA is code" idea.  At a more general level than whether specific genome structures match specific programming structures.
UPDATE #2: The late Scott Federhen of NCBI wrote an article comparing genomic replication with the lambda operator.
UPDATE #3: It might be the case that LISP was inspired by the genome.

I'm afraid you're going to get much the same reception here as you did in Biology. Please [edit] and point to specific features of the genome that, to you, are similar to specific features of LISP. What would the "func" be in the genome? What are the vars? The genome is usually thought of in terms of natural languages with words (codons) making sentences (genes) and punctuation (start and stop codons). And I fear you are vastly underestimating the complexity of the cellular machinery that reads the genome. Also, your first link is broken.

AKA Richard Dawkins without the proofs. You need to draw down a specific testable hypothesis, the days when rhetoric ruled evolutionary theory ended with sequencing. LISP has a lot of fans (Perl people like it), biological systems are far less tractable and much more diverse.

https://xkcd.com/793/

@terdon I explained exactly what the correspondence is and a simple prediction that results.

Well, operons are only relevant to bacterial genomes. What do you mean by "structural" genes? And yes, genes can be nested, but they can also overlap on opposite strands, are spliced, up and down regulated, all sorts of things. And then there are exceptions, edge cases, outliers... I still don't see what you're getting at here to be honest.

@terdon it looks like the genome uses a general structure that is the same as the basic LISP structure, so any theory regarding LISP type languages, and programming languages in general, is applicable to understanding the genome.  I'm asking if my understanding is correct, and if anyone has undertaken this sort of analysis before, i.e. analyzing the genome as if it is a programming language.  Let me know if this is still unclear and what more I can do to clarify my question.

I guess it may come down to my own ignorance (I am a biologist, not a computer scientist) but I have no idea what kind of theory you would be able to apply or how. My only contact with LISP is my .emacsrc file. But many of us here aren't engineers, so your question would benefit from a more specific explanation of what sort of tools/approaches you are thinking of. More to the point, I still feel you are not grasping just how complex and, most importantly, stochastic and chaotic system a genome is. I just don't see what benefit your approach would bring or how it would even be possible.

@terdon my basic idea is that if the LISP structure is consistent, then we could maybe look for other programming constructs built on the structure, such as recursion/looping, maybe templating or even object oriented programming.  Or just identify simple properties like balanced opening closing markers.  The mere fact it has a LISP structure would be fascinating if true.

@terdon also, more generally, I'm asking if anyone has taken this approach before, and demonstrated its effectiveness or found out there is some major problem so analyzing the genome as a programming language does not yield anything useful for a particular reason.  Do you happen to know either way what research has been done in this area?  Intuitively it seems a useful idea, because there is a wealth of knowledge regarding reverse engineering of programming languages, so we might be able to apply this knowledge to reverse engineer the genome.

Please [edit] your question and explain all this. Ideally, explain why you think programming languages are a good paradigm and better than natural languages. I still haven't seen any good evidence presented here and fear it might be a case of "when all you have is a hammer, everything looks like a nail". Granted, perhaps I'm the one with the hammer and not you, but that's why it would be really helpful if you could express your idea in terms of specific features of programming languages and how those would help you understand the genome better than we do already.

@terdon I have added the explanation from the comments to the question.  Hopefully it is clear enough now.

I think the name we use for "reverse engineering knowledge base to better understand the genome" is genome annotation.
But more seriously. I don't think there is any person who would possess sufficient knowledge of both to find practical analogies. Because we have plenty of ways how to analyze genome from different perspectives (e.g. genomes as subjects of mutations and selection over millions of years)m hence trying to fit some cs theories might be... difficult.

@KamilSJaron surely one more option might be helpful?  Why do you think this question is encountering so much resistance?  I've tried clarifying it a number of times, and in my opinion it is pretty clear what I am asking.

I'm voting to close this question as off-topic because it has nothing to do with bioinformatics

@MichaelG. seems like interpreting the genome as code is precisely bioinformatics.  Per the Wikipedia definition: As an interdisciplinary field of science, bioinformatics combines biology, computer science, information engineering, mathematics and statistics to analyze and interpret biological data.

@yters: bioinformatics is interpreting biological data with computers. The opposite is sometimes called "computational biology" which is using biological systems AS computers. Or analysing biological systems as if they were a computer. So DNA computing, for example.

@gilleain this is a bit different than using the systems as computers.  It is interpreting the systems as IT constructs, in this case a programming language.  Reading the Wikipedia articles on both topics, it still seems to me that bioinformatics is the best match, since the idea I propose relies on using computer science concepts to analyze the genome.  There is at least a surface visual similarity that seems quite strong, so I'm curious whether people have looked into this idea and what they have discovered.

I see what the metaphor has been inspired by, where I think it got turned around, and might have a better one.
I'd like to back up a bit and establish some common ground.  It isn't wrong to expect that some aspects of genetic code and software code are similar, and while biology is usually very specific, some aspects can be seen productively from the light of computer science ideas.  For example, mass-action molecular dynamics, which are often a good model for protein interactions, are Turing-complete.  Here's a fun example where people have written a programming language out of molecular dynamics:  https://arxiv.org/abs/1809.07430
It's also true that proteins have particular substructures leading to certain functions, such as using alpha-helices to react less with the outside world.  Learning about larger substructures, called protein domains, help identify the functions of families of proteins.
The problem is that the parts of the genome which appear nested like LISP cells actually aren't included in the protein.  The regulatory regions, promoters, inhibitors, introns, and so forth never become part of the key "computational" machinery.  This nesting isn't actively part of the computational activity proteins do perform.  The main quality of LISP is that the structure and function are perfectly self-similar throughout an expression, but the structure of the regulatory regions are neither replicated down into the protein nor are necessarily recapitulated on a larger scale.
Unfortunately, we don't expect the computational structure of proteins to be the same from one protein to another.  This is because "binds only with X" and "binds only with Y" must be implemented differently to successfully be distinct operations, and it is these binding constraints by which molecular dynamics implement "computational" operations, such as "not X" and "not Y", and so those are distinct as well.
I think a better way to think about the regular structures of the regulatory apparatus is more like a network protocol.  Along an otherwise undifferentiated binary channel, one has to identify when and where the data is delivered to an application.  Similarly, along a strand of bases, the machinery of translation has to have an indication of where the "application" of the protein starts and ends, and when that protein should be delivered to the cell.  We might compare some of the images of genetic structure that provided in the original question with layouts of network headers: http://slideplayer.com/slide/4798836/15/images/16/TCP/IP+Packet+Structures.jpg  We might compare bacterial versus multi-cellular regulation as having different delivery mechanisms and data structures for delivering the application information, for different functional consequences, in the same way TCP and UDP have different structures for different trade-offs.
This is all complicated somewhat by the fact that proteins do then interact with the regulation of genes, leading to interesting "computational" dynamics such as feed-forward regulation.  However, I'm not aware of any particular structure of the regulatory apparatus that necessarily correlates with these larger functional patterns.  The opposite is frequently biologically useful: the same structure of a protein helping to activate one gene might simultaneously inhibit another.
One book I've found very accessible and useful in understanding where computation is useful biologically is Dennis Bray's "Wetware: a computer in every living cell".  It gave me a clearer picture of how computation occurs in cells.  That should open up some introductions to systems biology which explains the biological-role of larger molecular "computational" patterns.

Great answer!  First answer of substance on either site!  I'll have to think it over, and if I find your reasoning persuasive, I'll accept this answer.

Thanks!  I'm glad you found the answer constructive.  After I found a useful perspective, I enjoyed thinking about the the question.

You say one problem with the idea is that the LISP like structure doesn't actually become part of the computational machinery.  I'm not sure this invalidates the connection I see.  When a LISP program is compiled down to the machine code format, the high level structure disappears.  This sounds like what you are saying with regard to proteins.  All the structural patterns disappear once the genome is compiled into a protein.  So, to my ears it seems you are saying LISP and the genome are the same in this regard, instead of different.  Perhaps you can clarify?  BTW, good point about protocols.

A high level comment, I think your network protocol idea and this LISP idea are consistent with a broader idea that the genome may have a number of different protocols that are followed, just like a modern distributed computing system.  And while they are understood holistically in reference to each other, each protocol can also be analyzed by itself.  Perhaps the same mode of thinking applies to the genetic code and regulation networks.  While they can be understood holistically, they perhaps can also be understood individually.

Right, so it is true that a kind of compilation does happen, where sequences of bases undergo a translation into the amino acids that make the protein.  The observation that I'm making is that the part of the gene in which you observed the LISP structure, is not included in that translation process.  Instead, those structures have been stripped away at that point.  This is similar to how when data is provided by a network protocol to the application, those headers are not delivered to the application.

The question is then whether or not we expect to see regular structures corresponding to computational structures in the part of the gene that is then translated, such as how we would see machine instructions for branching, adding, and these sorts of things.  Unfortunately, the answer is generally not, as the computations implied have to be specific for each different input.  If the answer were yes, the challenges of bioinformatics would be much easier!

It occurs to me that it might be helpful to know that promoters, inhibitors, introns, and so forth are occasionally called "non-coding regions" given that they aren't included in translation.  The following diagram might be helpful:  https://en.wikipedia.org/wiki/File:Gene_structure_prokaryote_2_annotated.svg

I see, just the genes inside the segment remain when translation takes place.  So, in that sense it is dissimilar to LISP, since the function in a LISP s-exp would be a part of the compilation process.  On the other hand, the non-coding section does seem to play a role in whether and how the genes are translated, if I'm understanding the answers so far.  And, thus this prokaryote genome structure still looks a lot like a LISP s-exp, since it can presumably also be nested: https://en.wikipedia.org/wiki/S-expression

As you've inferred, it is entirely right that these non-coding sections do play a role in regulation.  The way that regulation generally works, crudely speaking, is that whether transcription can go ahead or not depends upon how proteins bind to these regulatory sections, either guiding or blocking the transcription machinery.  The content of these regulatory sections must then be specific to the proteins that do work, so that once again the computational work that is done unfortunately must be specific to each family of proteins and not part of a more general computational mechanism.

I agree that there is a kind of nesting that can be represented as an S-expression.  Instead, it's my contention that the S-expression could be a misleading choice for representing that structure, as I propose that the nesting that you've observed is better explained as the nesting of protocols rather than the recursive nesting of context-free grammars.  By protocol nesting, I mean that after the data has been delivered from one protocol, and the metadata of that protocol is removed, the resulting data has the data and metadata of the next protocol.

Thanks again, this is a truly excellent series of answers, and you've clarified the issue greatly in my mind and given me a reference to further my understanding.  Great work!

Great, I'm glad to read that this has been helpful.  I have to hand it to you for not being dissuaded from pursuing a potential hypothesis for the structure you observed, while yet entertaining suggestions counter to your hypothesis, until you found an explanation that attempted to address all the issues.

So, first of all, is my high level observation correct? Does the genome, at least for bacteria, have a global LISP like structure?

No. I'm afraid this is what it boils down to. I don't see any real similarity to LISP or other programming languages at all. Perhaps it would be clearer if you could explain in detail why you think they're similar, but as it stands, no, there is no such structure apart from a very simplistic similarity which is sort of the same as saying that the genome is like a car because it carries genes and could therefore be studied by a car mechanic.

You don't find the correspondence list in my question clear enough?  I'm asking about the structure, which seems to be the same as a LISP program, or perhaps more generally an abstract syntax tree.  And yes, the structure itself is pretty simple, but still there seems to be an exact correspondence, more so than referring to a car. I am somewhat bemused this question is so hard for people to understand and answer.

@yters no,  I don't find it clear at all. As far as I can tell, you seem to be comparing functions and their arguments to operons and their genes, but I don't see any similarity at all there. That's why I posted so many comments asking you to get more specific. I have no idea why you feel that operons are more like functions than they are like cars or bowls or anything else that can contain items. You might want to look into gene regulatory networks though, those are closer to what you're getting at.

You could say there is a disanalogy in that the promoter/operator is not really a function, but controls whether the genes in the operon are "executed".  But there is still generally a notion of "execute this thing" and a nested tree structure.

@yters not really. That's not how the promoter works. Promoters tend to control how much it is executed, not whether. And they're just the tip of the iceberg in terms of regulation.

but is it correct the promoter/operator preamble somehow controls the manner in which the subsequent genes are 'executed'?

@yters it is only one of many things that affect gene expression, some of which are before the gene, some in it, some after it, some can even be quite distant from it and others aren't even part of the genome at all. As I said before, you are vastly underestimating the complexity of the genome.

I don't deny the extreme complexity of the greater context.  But, just as I can understand the C code of some UNIX utility without having to understand UNIX, perhaps we can analyze a genome by itself, to at least some degree.

You mention this as the base similarity you find:

operator = function
       structural genes = variables
       structural genes can be nested like cons cells

As other has pointed, this just considers the genome as a sequence of where only the linear order is important. It doesn't take into account other things that happen on the genome that are know to affect the cells like the chromatin structure, the duplication of certain regions, the transversion, the effect of moving elements, methylation, or simply that a single gene can have multiple functions (pleiotropism, moonlighting) while a programming function does not.
Even if we considered it as code, certainly by using tools like HMM we could also learn more about the structure of the genes we have. While certainly, we need more tools than just the sequence to discover what does each gene do.
There are a lot more differences than similarities on this comparison, it can be helpful for certain aspects but it fails for other aspects of cells (without even going into the virus discussion or the mitochondrial).

The only possible parallel was the "processor 'core' wars" in Avida and Tierra. Tierra was initiated by Tom Ray, around 1990. Here they had a 'population' of compiled binaries were randomly mutated (yep the binary code was mutated) and reproduction was dependent on how much of the processor could be consumed. The simulation worked and an in silico ecosystem evolved, but was heavily dependent on the preconditions of the simulation and proved no more reliabe at understanding evolution than differential equestions. So it was a nice idea, and very cool it worked, but it came to zero.
That aside the best example was the concept of the "selfish gene", but has very little to do with programming languages IMO and everything to do with understanding kinselection.

Not sure what question you are answering here.  I am not asking about evolution simulations.

To say I want "genome understanding" without evolution, is like saying I want a watch without hands. Can't have one without the other.

I would say it is more like I can understand the C source code of a UNIX utility without having to know the history of UNIX or C's development.

Believe it or not, there are people who know more than you do about programming, and who also know a hell of a lot more than you do about the genome.  This is not an avenue those people think is worth pursuing.
DNA is a chemical, and interacts with other chemicals based on its shape.  There is no equivalent to that in a human designed programming language.

This seems unnecessarily combative. The OP never claimed to be an expert!

Which is why OPs first step should have been to read what experts in both subjects have been doing in the field he's interested in, instead of trying to impress a bunch of volunteers with programming jargon.

Wow. Please take a step back and stop assuming malice. The OP did do some research and did try and see if any work had been done. Nobody's out to impress anybody here, so please stop being so aggressive!


text:
Raw Data frame manipulation in python
Using python 3 I need to process qPCR sequencing raw data outputs by searching for the first occurrence of a user defined string and then making a new data frame using all lines after that string. I am trying to find solutions in the pandas doc but so far unsuccessful.
This is a raw output .csv file that I need to process. (couldn't paste complete csv as exceeds character limit, this is lines 40-50 and am hoping this text is useful?). I need to tell pandas to create a new data frame that 1. starts at the line containg the first occurance of str("Sample Name") with that line as header and containing all lines following. And then 2., only including columns ("Sample Name"), ("Target Name"), ("CT").
Could someone please help me so that I can use python to help me analyze biological data?
Many thanks,
Luke
40,Quantification Cycle Method,Ct,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

41,Signal Smoothing On,true,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

42,Stage where Melt Analysis is performed,Stage3,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

43,Stage/ Cycle where Ct Analysis is performed,"Stage2, Step2",,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

44,User Name,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

45,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

46,Well,Well Position,Omit,Sample Name,Target Name,Task,Reporter,Quencher,Quantity,Quantity Mean,SE,RQ,RQ Min,RQ Max,CT,Ct Mean,Ct SD,Delta Ct,Delta Ct Mean,Delta Ct SD,Delta Ct SE,Delta Delta Ct,Automatic Ct Threshold,Ct Threshold,Automatic Baseline,Baseline Start,Baseline End,Amp Status,Comments,Cq Conf,CQCONF,HIGHSD,OUTLIERRG,Tm1,Tm2,Tm3,Tm4

47,1,A1,False,WT1,AtTubulin,UNKNOWN,SYBR,None,,,,,,,23.357698440551758,23.4766845703125,0.5336655378341675,,,,,,True,20959.612776965325,True,3,17,Amp,,0.9588544573203085,N,Y,N,81.40960693359375,,,

48,2,A2,False,WT1,AtTubulin,UNKNOWN,SYBR,None,,,,,,,24.05980110168457,23.4766845703125,0.5336655378341675,,,,,,True,20959.612776965325,True,3,15,Amp,,0.9592687354496955,N,Y,N,81.40960693359375,,,

49,3,A3,False,WT1,AtTubulin,UNKNOWN,SYBR,None,,,,,,,23.012556076049805,23.4766845703125,0.5336655378341675,,,,,,True,20959.612776965325,True,3,16,Amp,,0.9592714462250367,N,Y,N,81.40960693359375,,,

50,4,A4,False,fla11fla12-1,AtTubulin,UNKNOWN,SYBR,None,,,,,,,23.803699493408203,24.419523239135742,0.5669151544570923,,,,,,True,20959.612776965325,True,3,17,Amp,,0.9671570584141241,N,Y,N,81.40960693359375,,,

This is the code that I have so far:
import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

data = pd.read_excel ("2019-02-27_161601 AtWAKL8 different version expressions.xls", sheet_name='Results').fillna(0)

data.to_csv('df1' + '.csv', index=True)

df1 = pd.read_csv ("df1.csv")

Please [edit] your question and add i) an example of your input data, ii) an explanation of the data frame you want to build from that data and iii) the code you have so far so we don't reinvent the wheel.

Don't use grep using subprocess to filter dataframes. That's crazy :-D
Follow a pandas tutorial: the filtering that you want to do on your dataframe can be done with pandas commands, no need to use grep. If you are using python, avoid mixing your code with shell commands.

Please don't give us photographs of your data! We need your actual data so we can test any solutions we come up with. Just save the excel file as csv or tsv or something and post the first few lines here. And paste the errors as well. You can use the formatting tools to make all this look like code.

What is that header about in your csv file? I suppose that a csv file should be just a comma-separated table. What is the first column? line numbers? Why do they start at 32? I understand that starting with python might be difficult, have you tried some tutorials?

Thanks for all resposes. @terdon that formatting link helped,

@Kamil S Jaron I have tried to explain my csv example and what I need to acheive a bit better, hope that's helpful?

data = pd.read_excel ("2019-02-27_161601 AtWAKL8 different version expressions.xls", sheet_name='Results').fillna(0)

data.to_csv('df1' + '.csv', index=True)

df1 = pd.read_csv ("df1.csv")

This doesn't make much sense ... you read in Excel to a dataframe, you then read out to csv - the dataframe doesn't have an index, but you include an index and then read it in again. This could do more harm than good.
Anyway your supervisor is right it is very easy to parse this data using pandas, although there will be more expertise doing this via R.
I don't want to peer through the gloom of a CSV output, but I would like to highlight how easy it is to parse this stuff into a dataframe. The dataframe 'data' is fine.

Parse You can simply remove the top rows of bumf data = data.iloc[6:] (I might have miss counted, plus see below),
Headers You will probably have to remove the "headers" row to a list and
then introduce back onto the dataframe using data.columns =
['header1', 'header2' ...] ... rather this is the list (array in Perl) which you can automatically generate using a magic to_list type pandas command which I forget, its not list - its cuter
Columns You next might use drop to remove the columns you don't want (that might be alot of columns) or else make a subset of
the dataframe into a new dataframe, you can also use iloc for
columns too.
Index You then assign A1, A2, ... etc  as an index data.index = ['A1',
'A2', ]...etc..  you can reshape column 3 to automate this
to_numeric You might have to assign the CT values as floats
    or better to_numeric
Extrapolation You would then define a subroutine (or whatever they call them in
Python) so perform the extrapolation and use the apply command on
your dataframe and assign it to a new column, where the input is your CT values.
Output Finally, issue your to_csv command to_pickle is also a good idea because you have a df that is worth importing at a future date.
Plotting this is easy using data.plot which appears to be some sort of inbuilt matplotlib, keep in mind seaborn DONT try using R's ggplot2 within Python using rpy2 .... its a NIGHTMARE to get working and doesn't support the latest ggplot2 (!!!). R
s ggplot2 is stronger  than seaborn or matplotlib (providing you are working within R), but these two are good for quick graphs.

If I had time I'd write the code for you, but I'm afraid I'll have to leave it there. I'm still in debt with a Perl efetch script I promised a previous OP.
In summary, Python pandas is fabulous stuff for this purpose, I love it. Doing this via Perl and "subprocessing" inside R... its doable, well in this particular case I suspect its just as easy, but generally its alot more code and often fiddly.

Really appreciate your help! Thanks. When I get a chance I will post my updated code.

authors: Олеся Галанова


text:
Subset a protein fasta based on sequence length?
I have a strange protein fasta file. There are several entries for the same gene and I need to extract the entry with the longest width.
>lcl|NW_017095468.1_prot_XP_017786561.1_8 [gene=LOC108569215] [db_xref=GeneID:108569215] [protein=NADH dehydrogenase [ubiquinone] 1 a
lpha subcomplex subunit 9, mitochondrial] [protein_id=XP_017786561.1] [location=complement(join(88030..88986,89044..89243,89412..8944
8))] [gbkey=CDS]
MAAVIFTGAQLLKQQSGLVGIAYIRVNNYSSDAKVYNLASLKRGTGGRSSFNGIVATVFGSTGFIGRYVANKLGKIGTQL
ILPYRGNNYETMRLKLCGDLGQVLFQPYYLKDLESIDKALKYSNVVINLVGRDWETKNFSFQDVHVKGARDLARAAKKAG
VEKFIHLSALNCDGPNEAIFSRQGSKFLSSKWEGEQAVLEEFPEATIFRPSDVYGQEDRFLRYYAHNWRRQGQYMPMWKN
GEATIKQPVHVSDLAAGIVAAIKDPEAAGKVYQAVGPRRYQLNELVDWFYRVMRKDSEWGYKRYDMKYDMFFKLKVSLTQ
KFSPAFPVGNLHWERLEREFVTDRVNQAIPTLEDLGVNLRRMEDQVPWELKPFTYGLYFGTDSEEPVVEPKPPKYVS
>lcl|NW_017095468.1_prot_XP_017770525.1_9 [gene=LOC108557919] [db_xref=GeneID:108557919] [protein=spondin-1 isoform X1] [protein_id=XP_017770525.1] [location=join(95471..95648,106560..106774,106834..107025,107093..107299,107562..107816,107872..108474,108526..108799,108885..109058,109127..109237,109573..109812,109882..110114)] [gbkey=CDS]
MRLKVAFLWLVSSISWIGEALRCDRTPEGTFSPRTRADGRFVIEVSGNPDTYVPGEQYNIFLRSNGEYQAKNKFKDFLLL
VEHEPSEKILGEVHNPSVGTLQLLGDMLMKFSEKCRNAVMQTNSLPKSEVQVLWVAPPSGSGCVAIRATVVESKEFWYTD
DGPLSKILCEEVQENEDTQPNMLRQCCACDEAKYEVTFEGLWSRNTHPKDFPSNGWLTRFSDIIGASHTFDYTFWNYGEI
ASNGLRQLAENGNTRMLESELKAKSEHIRTIIKARGISYPNITGKTFAVFRVDKRHHLMSLVSMIDPSPDWIVGVSGLEL
CLRNCSWVESRVLNLYPWDAGTDDGPTYISANQPSMPPHPIRRIKSNSPNDPRSPFYDPTGTEMKPLARLYLSRQRLYEK
NCVAQVDVSEEDGGVEGDKCEMEEWSEWSKCTVTCGRGFKYKQRAYKNPASNFVCNKPLTKRASCVAILEHCSNQQRPQE
ADPSCSLTGWGNWSSCTAPCGPGWKTRSRRYKNRKAAKRCAAGNENPEPLEQNLECMERECGPSDRRPLQESKECEGRAW
SNWSPCSSTCGKGIKVRRRMAYRSLWGRSPARYNRGLFDTEDTSRDDDDDGSDEDPCMNLDEKVECINDDVPVCEDTVDN
SAVVCGFPRDEGGCMSNVDRWYFDVIKGNCDIFSYSGCQGNKNNFKTLERCENVCDSYKKELLANRTAYKRQLGVTVSGV
LSYNLHHMQNDDADNCVPGSQTRQDQDKKIIQEPIGEVVDCQMSEWTNWSGCNATCGRGFSTKHRFIRVHPSNGGKRCPQ
KILKKRKCKIPCPGDYTKRDPMLPTWGTANSLEHVQIDCVMTGWSAWSPCSRSCGPNAVQQRTRGILLPPSGRGEPCLHR
TEERPCSLLACPE
>lcl|NW_017095468.1_prot_XP_017771299.1_10 [gene=LOC108557919] [db_xref=GeneID:108557919] [protein=spondin-1 isoform X1] [protein_id=XP_017771299.1] [location=join(95471..95648,106560..106774,106834..107025,107093..107299,107562..107816,107872..108474,108526..108799,108885..109058,109127..109237,109573..109812,109882..110114)] [gbkey=CDS]
MRLKVAFLWLVSSISWIGEALRCDRTPEGTFSPRTRADGRFVIEVSGNPDTYVPGEQYNIFLRSNGEYQAKNKFKDFLLL
VEHEPSEKILGEVHNPSVGTLQLLGDMLMKFSEKCRNAVMQTNSLPKSEVQVLWVAPPSGSGCVAIRATVVESKEFWYTD
DGPLSKILCEEVQENEDTQPNMLRQCCACDEAKYEVTFEGLWSRNTHPKDFPSNGWLTRFSDIIGASHTFDYTFWNYGEI
ASNGLRQLAENGNTRMLESELKAKSEHIRTIIKARGISYPNITGKTFAVFRVDKRHHLMSLVSMIDPSPDWIVGVSGLEL
CLRNCSWVESRVLNLYPWDAGTDDGPTYISANQPSMPPHPIRRIKSNSPNDPRSPFYDPTGTEMKPLARLYLSRQRLYEK
NCVAQVDVSEEDGGVEGDKCEMEEWSEWSKCTVTCGRGFKYKQRAYKNPASNFVCNKPLTKRASCVAILEHCSNQQRPQE
ADPSCSLTGWGNWSSCTAPCGPGWKTRSRRYKNRKAAKRCAAGNENPEPLEQNLECMERECGPSDRRPLQESKECEGRAW
SNWSPCSSTCGKGIKVRRRMAYRSLWGRSPARYNRGLFDTEDTSRDDDDDGSDEDPCMNLDEKVECINDDVPVCEDTVDN
SAVVCGFPRDEGGCMSNVDRWYFDVIKGNCDIFSYSGCQGNKNNFKTLERCENVCDSYKKELLANRTAYKRQLGVTVSGV
LSYNLHHMQNDDADNCVPGSQTRQDQDKKIIQEPIGEVVDCQMSEWTNWSGCNATCGRGFSTKHRFIRVHPSNGGKRCPQ
KILKKRKCKIPCPGDYTKRDPMLPTWGTANSLEHVQIDCVMTGWSAWSPCSRSCGPNAVQQRTRGILLPPSGRGEPCLHR
TEERPCSLLACPE
>lcl|NW_017095468.1_prot_XP_017772069.1_11 [gene=LOC108557919] [db_xref=GeneID:108557919] [protein=spondin-1 isoform X2] [protein_id=XP_017772069.1] [location=join(95471..95648,106560..106774,106834..107025,107093..107299,107562..107816,107872..108474,108526..108799,109546..109812,109882..110114)] [gbkey=CDS]
MRLKVAFLWLVSSISWIGEALRCDRTPEGTFSPRTRADGRFVIEVSGNPDTYVPGEQYNIFLRSNGEYQAKNKFKDFLLL
VEHEPSEKILGEVHNPSVGTLQLLGDMLMKFSEKCRNAVMQTNSLPKSEVQVLWVAPPSGSGCVAIRATVVESKEFWYTD
DGPLSKILCEEVQENEDTQPNMLRQCCACDEAKYEVTFEGLWSRNTHPKDFPSNGWLTRFSDIIGASHTFDYTFWNYGEI
ASNGLRQLAENGNTRMLESELKAKSEHIRTIIKARGISYPNITGKTFAVFRVDKRHHLMSLVSMIDPSPDWIVGVSGLEL
CLRNCSWVESRVLNLYPWDAGTDDGPTYISANQPSMPPHPIRRIKSNSPNDPRSPFYDPTGTEMKPLARLYLSRQRLYEK
NCVAQVDVSEEDGGVEGDKCEMEEWSEWSKCTVTCGRGFKYKQRAYKNPASNFVCNKPLTKRASCVAILEHCSNQQRPQE
ADPSCSLTGWGNWSSCTAPCGPGWKTRSRRYKNRKAAKRCAAGNENPEPLEQNLECMERECGPSDRRPLQESKECEGRAW
SNWSPCSSTCGKGIKVRRRMAYRSLWGRSPARYNRGLFDTEDTSRDDDDDGSDEDPCMNLDEKVECINDDVPVCEDTVDN
SVTENYFGRIVPGSQTRQDQDKKIIQEPIGEVVDCQMSEWTNWSGCNATCGRGFSTKHRFIRVHPSNGGKRCPQKILKKR
KCKIPCPGDYTKRDPMLPTWGTANSLEHVQIDCVMTGWSAWSPCSRSCGPNAVQQRTRGILLPPSGRGEPCLHRTEERPC
SLLACPE
>lcl|NW_017095468.1_prot_XP_017772815.1_12 [gene=LOC108557919] [db_xref=GeneID:108557919] [protein=spondin-1 isoform X3] [protein_id=XP_017772815.1] [location=join(95471..95648,106560..106774,106834..107025,107093..107299,107562..107816,107872..108474,108526..108799,109573..109812,109882..110114)] [gbkey=CDS]
MRLKVAFLWLVSSISWIGEALRCDRTPEGTFSPRTRADGRFVIEVSGNPDTYVPGEQYNIFLRSNGEYQAKNKFKDFLLL
VEHEPSEKILGEVHNPSVGTLQLLGDMLMKFSEKCRNAVMQTNSLPKSEVQVLWVAPPSGSGCVAIRATVVESKEFWYTD

I've tried using BioStrings but I'm unfamiliar with how to index this type of file. Any help appreciated.

Can you give us an example that includes a duplicated gene? How is the "gene" defined? Should we look for the RefSeq protein ID (XP_123), for the gene=LOC string? All of the above? Will the entire header line be identical or only parts of it?

I've added more of the file to my question. I think the unique/isoform bits are the XP_ (which has information about isoforms) ... NW_ bits are the genes.

NW_017095468 is a whole genome shotgun sequence, it will have many genes. The XP bits are RefSeq protein accessions. These can be different for different protein products of the same gene. It looks like the solutions will need to look at either the gene= or the GeneID: values.

Perhaps more important is for R-MASHup to explain what is "strange" about these sequences and why. It is well known that a gene may encode multiple proteins, so the fact that there are multiple entries per gene doesn't surprise me.

Yes. And choosing the longest is kind of arbitrary, but all this depends on what you're doing. For instance, I once needed to do something very similar because I was looking at various features of certain genes and I wanted to see if protein product length correlated with some other characteristic, so I chose the longest isoform of each gene's products. But, @R-MASHup, if you think the longest will be somehow special (most highly expressed or more representative or anything like that), you will be surprised.

Well it's a nonmodel organism, trying to find orthologs, and so need to reduce the number of entries somehow (esp as it's computationally so slow) ... how would you pick?

For orthologs, I'd probably not pick since it's impossible to predict which protein will have ended up shared with the other organism. I'd just run a tblastn or better exonerate on using all proteins as queries. But yeah, picking the longest might make sense if you don't have a decent machine. Also, if my answer isn't working for you, you probably have ID lines that are different from what you've shown here. My solution just takes the 2nd field of the line ([gene=LOCNNNN] in your example), if that changes, it will not work.

Using the FastaToTbl and TblToFasta scripts I have posted before, you can do:
FastaToTbl file.fa |
    awk '{
            len=length($NF);
            if(len > l[$2]){
                l[$2]=len;
                seq[$2]=$0;
            }
          }
          END{
            for(s in seq){
                print seq[s]
            }
        }' | TblToFasta

Running this on your example gives:
$ FastaToTbl file.fa | awk '{len=length($NF); if(len > l[$2]){l[$2]=len; seq[$2]=$0; }} END{for(s in seq){print seq[s]}}' | TblToFasta
>lcl|NW_017095468.1_prot_XP_017770525.1_9 [gene=LOC108557919] [db_xref=GeneID:108557919] [protein=spondin-1 isoform X1] [protein_id=XP_017770525.1] [location=join(95471..95648,106560..106774,106834..107025,107093..107299,107562..107816,107872..108474,108526..108799,108885..109058,109127..109237,109573..109812,109882..110114)] [gbkey=CDS]
MRLKVAFLWLVSSISWIGEALRCDRTPEGTFSPRTRADGRFVIEVSGNPDTYVPGEQYNI
FLRSNGEYQAKNKFKDFLLLVEHEPSEKILGEVHNPSVGTLQLLGDMLMKFSEKCRNAVM
QTNSLPKSEVQVLWVAPPSGSGCVAIRATVVESKEFWYTDDGPLSKILCEEVQENEDTQP
NMLRQCCACDEAKYEVTFEGLWSRNTHPKDFPSNGWLTRFSDIIGASHTFDYTFWNYGEI
ASNGLRQLAENGNTRMLESELKAKSEHIRTIIKARGISYPNITGKTFAVFRVDKRHHLMS
LVSMIDPSPDWIVGVSGLELCLRNCSWVESRVLNLYPWDAGTDDGPTYISANQPSMPPHP
IRRIKSNSPNDPRSPFYDPTGTEMKPLARLYLSRQRLYEKNCVAQVDVSEEDGGVEGDKC
EMEEWSEWSKCTVTCGRGFKYKQRAYKNPASNFVCNKPLTKRASCVAILEHCSNQQRPQE
ADPSCSLTGWGNWSSCTAPCGPGWKTRSRRYKNRKAAKRCAAGNENPEPLEQNLECMERE
CGPSDRRPLQESKECEGRAWSNWSPCSSTCGKGIKVRRRMAYRSLWGRSPARYNRGLFDT
EDTSRDDDDDGSDEDPCMNLDEKVECINDDVPVCEDTVDNSAVVCGFPRDEGGCMSNVDR
WYFDVIKGNCDIFSYSGCQGNKNNFKTLERCENVCDSYKKELLANRTAYKRQLGVTVSGV
LSYNLHHMQNDDADNCVPGSQTRQDQDKKIIQEPIGEVVDCQMSEWTNWSGCNATCGRGF
STKHRFIRVHPSNGGKRCPQKILKKRKCKIPCPGDYTKRDPMLPTWGTANSLEHVQIDCV
MTGWSAWSPCSRSCGPNAVQQRTRGILLPPSGRGEPCLHRTEERPCSLLACPE
>lcl|NW_017095468.1_prot_XP_017786561.1_8 [gene=LOC108569215] [db_xref=GeneID:108569215] [protein=NADH dehydrogenase [ubiquinone] 1 alpha subcomplex subunit 9, mitochondrial] [protein_id=XP_017786561.1] [location=complement(join(88030..88986,89044..89243,89412..89448))] [gbkey=CDS]
MAAVIFTGAQLLKQQSGLVGIAYIRVNNYSSDAKVYNLASLKRGTGGRSSFNGIVATVFG
STGFIGRYVANKLGKIGTQLILPYRGNNYETMRLKLCGDLGQVLFQPYYLKDLESIDKAL
KYSNVVINLVGRDWETKNFSFQDVHVKGARDLARAAKKAGVEKFIHLSALNCDGPNEAIF
SRQGSKFLSSKWEGEQAVLEEFPEATIFRPSDVYGQEDRFLRYYAHNWRRQGQYMPMWKN
GEATIKQPVHVSDLAAGIVAAIKDPEAAGKVYQAVGPRRYQLNELVDWFYRVMRKDSEWG
YKRYDMKYDMFFKLKVSLTQKFSPAFPVGNLHWERLEREFVTDRVNQAIPTLEDLGVNLR
RMEDQVPWELKPFTYGLYFGTDSEEPVVEPKPPKYVS


text:
Cluster is split in 2-3 locations on tsne plot - Suerat
I am running a single cell dataset (count data - exon) through Seurat. After running tsne I see a cluster (13) split in 3 different locations on the plot. Here are the commands I am running:
library('Seurat')
## rows (genes) that have 99% zeros were filtered out,
## and then columns (cells) that have 99% zeros were filtered out
isp.data <- read.table(file = "genes.cells99.percent.filtered", sep =
            ",", row.names = 1, header=TRUE)
isp <- CreateSeuratObject(raw.data = isp.data, project =
       "brain_regions")
isp <- NormalizeData(object = isp, normalization.method =
       "LogNormalize",scale.factor = 10000)
isp <- FindVariableGenes(object = isp, mean.function = ExpMean,
       dispersion.function = LogVMR,do.plot = FALSE)
hv.genes <- head(rownames(isp@hvg.info), 1000)
isp <- ScaleData(object = isp, genes.use = hv.genes, num.cores = 3)
isp <- RunPCA(object = isp, pc.genes = hv.genes, pcs.compute = 100,
       do.print = TRUE, pcs.print = 1:5, genes.print = 5)
PCElbowPlot(object = isp, num.pc = 100)
isp <- JackStraw(object = isp, num.replicate = 100, display.progress =
                 FALSE, num.pc = 30)
JackStrawPlot(object = isp, PCs = 1:30)
isp <- FindClusters(object = isp, reduction.type = "pca", dims.use =
        1:26, resolution= 0.75, print.output = 0, save.SNN = TRUE,
        force.recalc=TRUE)
isp <- RunTSNE(isp, dims.use = 1:26, do.fast = T)
TSNEPlot(isp, do.label = T, pt.size = 0.5)

Here is the tsne plot:

I have tried varying parameters for resolution in FindClusters (from 0.5-3.5), but that just increases the clusters found. Further I tried to play around with parameters min.dist (0.75,1), n_neighbors (10 to 30) and perplexity(5 to 50) but I still get at least 1-2 clusters split in the tsne plot.
Another thing I have tried is to reduce the dimension (dims.use=10) on which tsne is ran but the clusters seem to be too close together, like so:

Note: Dimensionality of dataset (after filtering) - 23824 genes across 8967 samples.
My question being why does such behavior happen? Is my filtering strategy not correct? Have I not normalized/scaled data correctly?

Your cluster labels come from graph clustering implemented in the FindClusters() function. The resulting clusters are then visualised with a 2D tSNE plot (via RunTSNE() and TSNEPlot()). So, your cluster 13 is not split into three sub-clusters but cells within cluster 13 look somewhat distant from each other on the tSNE plot. Having mentioned distances on a tSNE plot, please see this SE Data Science entry Can closer points be considered more similar in T-SNE visualization? and especially the comment "What I mean is that you can't just use distance in the lower space as a similarity criterion". In your case "the lower space" would be your tSNE plot.
As far as I can see from your code, you are following the widely-adopted "Seurat workflow" and all seem fine including your TSNE plot output. I could not spot the QC step though: apart from the one at object initiation, one usually applies filters for number of genes and % mitochondrial gene contribution (using FilterCells() in Seurat v2 and subset() in Seurat v3).
It would be a good idea to try UMAP instead of tSNE, please see Evaluation of UMAP as an alternative to t-SNE for single-cell data (Becht et al, 2019). Seurat v3 vignettes are now using UMAP and not tSNE any more.
Lastly, please see How to Use t-SNE Effectively to learn about points to take into account when using/interpreting tSNE plots. You will see that you would get surprisingly dissimilar tSNE plots from the same underlying data with different parameters (some are hardcoded in the RunTSNE() function).

Yes, I was missing step where I needed to filter for mitochondria and number of genes.

Great and very comprehensive answer!

Hi user2998764, I was just curious, do the tSNE plot looks "better" after the QC for % mito and genes?