text:
question to pick values based on conditions using a loop
I have a dataset with following information. Column C1_1 is chromosome number, C1_2 is SNP position and C3 is the p-value. I want to pick the most significant association within a genomic region of 3000 bp:
C1  C2  C3
L01_005000  L002g034    1.5e-12
L01_003000  L001g045    2.3e-6
L02_145000  L004g034    1.1e-2
L02_155000  L002g050    1.1e-4
L02_148000  L002g001    1.1e-3

I want to split first column into two as:
C1_1    C1_2    C2  C3
L01 005000  L002g034    1.5e-12
L01 003000  L001g045    2.3e-6
L02 145000  L004g034    1.1e-2
L02 155000  L002g050    1.1e-4
L02 148000  L002g001    1.1e-3

Then I want to sort subsets such as L01 and L02 separately based on C1_2 column:
C1_1    C1_2    C2  C3
L01 003000  L001g045    2.3e-6
L01 005000  L002g034    1.5e-12
L02 145000  L004g034    1.1e-2
L02 148000  L002g001    1.1e-3
L02 155000  L002g050    1.1e-4

Then I want to pick one value from each subset such as for L02 – I want to pick the most significant value based on C3 column under the difference of 3000. For example, L02-148000 is chosen as shown below as its C3 value is lower than L02-145000, and then the next value L02-155000 is chosen as the difference is more than 3000 than the previous value of L02-148000:
C1_1    C1_2    C2  C3
L01 005000  L002g034    1.5e-12
L02 148000  L002g001    1.1e-3
L02 155000  L002g050    1.1e-4

I have a long dataset with 500,000 rows of values and the subset goes from L01-L24. I think it will need a loop to process these many values to get the final values.
This is just an example, I have shown one for L01 because the difference in L01 203000 and L01 205000 is less than 3000. I have shown two values for L02 because the difference between L02 145000 and L02 148000 is not more than 3000. Whereas the difference between L02 155000 and L02 148000 is more than 3000.

Is your result C1_2 value for L01 supposed to be 005000 or 205000?

I think that is a typo @RamRS L01 should be 205000 (lower table), either that or its not possible to solve.

If @RamRS answers solves points 1 to 3, I'd recommend accepting this as the answer and then asking a further question on the 3000 threshold. In Python you'd need to construct a new column based the differences between C1_2 .. well thats how I'd do it anyway.

This question is deceptively challenging, but here is a potential solution:
library(tidyverse)

df <- read.table(text = "C1  C2  C3
L01_005000  L002g034    1.5e-12
L01_003000  L001g045    2.3e-6
L02_145000  L004g034    1.1e-2
L02_155000  L002g050    1.1e-4
L02_148000  L002g001    1.1e-3",
header = TRUE)

df
#>           C1       C2      C3
#> 1 L01_005000 L002g034 1.5e-12
#> 2 L01_003000 L001g045 2.3e-06
#> 3 L02_145000 L004g034 1.1e-02
#> 4 L02_155000 L002g050 1.1e-04
#> 5 L02_148000 L002g001 1.1e-03

df %>%
  separate(C1, into = c("C1_1", "C1_2")) %>%
  mutate(C1_2 = as.numeric(C1_2)) %>%
  arrange(C1_1, C1_2) %>%
  group_by(C1_1) %>%
  mutate(subgroup = cumsum(C1_2 > lag(C1_2, default = first(C1_2)) + 3000)) %>%
  group_by(C1_1, subgroup) %>%
  slice_min(C3) %>%
  ungroup() %>%
  select(-subgroup)
#> # A tibble: 3 × 4
#>   C1_1    C1_2 C2            C3
#>   <chr>  <dbl> <chr>      <dbl>
#> 1 L01     5000 L002g034 1.5e-12
#> 2 L02   148000 L002g001 1.1e- 3
#> 3 L02   155000 L002g050 1.1e- 4

## additional example data
df <- read.table(text = "C1  C2  C3
L01_005000  L002g034    1.5e-12
L01_003000  L001g045    2.3e-6
L02_145000  L004g034    1.1e-2
L02_155000  L002g050    1.1e-4
L02_158000  L004g034    1.1e-2
L02_159000  L002g050    1.1e-4
L03_165000  L002g001    1.1e-3",
header = TRUE)

df
#>           C1       C2      C3
#> 1 L01_005000 L002g034 1.5e-12
#> 2 L01_003000 L001g045 2.3e-06
#> 3 L02_145000 L004g034 1.1e-02
#> 4 L02_155000 L002g050 1.1e-04
#> 5 L02_158000 L004g034 1.1e-02
#> 6 L02_159000 L002g050 1.1e-04
#> 7 L03_165000 L002g001 1.1e-03

df %>%
  separate(C1, into = c("C1_1", "C1_2")) %>%
  mutate(C1_2 = as.numeric(C1_2)) %>%
  arrange(C1_1, C1_2) %>%
  group_by(C1_1) %>%
  mutate(subgroup = cumsum(C1_2 > lag(C1_2, default = first(C1_2)) + 3000)) %>%
  group_by(C1_1, subgroup) %>%
  slice_min(C3) %>%
  ungroup() %>%
  select(-subgroup)
#> # A tibble: 5 × 4
#>   C1_1    C1_2 C2            C3
#>   <chr>  <dbl> <chr>      <dbl>
#> 1 L01     5000 L002g034 1.5e-12
#> 2 L02   145000 L004g034 1.1e- 2
#> 3 L02   155000 L002g050 1.1e- 4
#> 4 L02   159000 L002g050 1.1e- 4
#> 5 L03   165000 L002g001 1.1e- 3

Created on 2023-02-07 with reprex v2.0.2
Edit
If you have two SNPs with the same C3 value (i.e. the top two rows in the example dataset below) you can select one of them using slice_min(C3, with_ties = FALSE), e.g. the code below will select the 'first' entry when two C3 values are the same within the same group/subgroup:
library(tidyverse)

df <- read.table(text = "C1  C2  C3
L01_005000  L002g034    1.5e-12
L01_003000  L001g045    1.5e-12
L02_145000  L004g034    1.1e-2
L02_155000  L002g050    1.1e-4
L02_158000  L004g034    1.1e-2
L02_159000  L002g050    1.1e-4
L03_165000  L002g001    1.1e-3",
header = TRUE)

df
#>           C1       C2      C3
#> 1 L01_005000 L002g034 1.5e-12
#> 2 L01_003000 L001g045 1.5e-12
#> 3 L02_145000 L004g034 1.1e-02
#> 4 L02_155000 L002g050 1.1e-04
#> 5 L02_158000 L004g034 1.1e-02
#> 6 L02_159000 L002g050 1.1e-04
#> 7 L03_165000 L002g001 1.1e-03

df %>%
  separate(C1, into = c("C1_1", "C1_2")) %>%
  mutate(C1_2 = as.numeric(C1_2)) %>%
  arrange(C1_1, C1_2) %>%
  group_by(C1_1) %>%
  mutate(subgroup = cumsum(C1_2 > lag(C1_2, default = first(C1_2)) + 3000)) %>%
  group_by(C1_1, subgroup) %>%
  slice_min(C3, with_ties = FALSE) %>%
  ungroup() %>%
  select(-subgroup)
#> # A tibble: 4 × 4
#>   C1_1    C1_2 C2            C3
#>   <chr>  <dbl> <chr>      <dbl>
#> 1 L01     3000 L001g045 1.5e-12
#> 2 L02   145000 L004g034 1.1e- 2
#> 3 L02   155000 L002g050 1.1e- 4
#> 4 L03   165000 L002g001 1.1e- 3

Created on 2023-02-08 with reprex v2.0.2

Wonderful use of cumsum + lag. Very nifty solution!

Thank you @jared_mamrot! I tried using your code. It gives the correct result with my data. Except one thing: if there are two SNPs with different positions in C1_2 but same C3 values within the 3000 window, it lists all those associations. I have not thought about this situation earlier.

Sounds like an easy fix @user1567654 - instead of slice_min(C3) try slice_min(C3, with_ties = FALSE) - that should do the trick, but if not, leave a comment and I'll look into it :)

EDIT: Ignore this answer, it does not address one of your requirements - that of the 3000-long bin. You're going to need some partitioning operations which are not straightforward.
@jared_marmot has cracked it in their solution.
This should be easy to do with dplyr/tidyr. Use
## Code not tested, might need tweaking
library(tidyr)
library(dplyr)
df_split <- df %>% separate(C1, sep = "_", into = c("C1_1","C1_2"))
df_aggr <- df_split %>% select(C1_1, C3) %>% group_by(C1_1) %>% summarise(C3 = min(C3)) %>% inner_join(df_split)

There are two problems that make your 3000 window complicated:

The start and end are not fixed: it's not 1 to n in steps of 3000; that'd be a straightforward bin.
We cannot even go with the smallest coordinate as the start, as you may have data points at 1500, 2500, 7200, 7600. If those data points exist, a set of 3000-long windows starting at 1500 would be 1500-4500, 4500-7500, etc. making your data go in 3 bins whereas your data points belong in 2 bins. You need windows with independent start coordinates. I've never encountered this requirement in my career even as a puzzle so I'm starting to think this might be an X-Y problem that needs to be addressed at a different level.

Oh wait you have that 3000 condition. That's not easy to implement. You'll need to create a new column with a unique value per 3000-bin that can be used to group bins together, which is not simple as the 3000's bounds are not fixed.

I think this is basically the answer and the 3000 threshold is a separate question.


text:
Generating sequences from HMM
I encountered this question in a textbook (source: https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.721.2540&rep=rep1&type=pdf) and was trying to workout the solution for it. Was wondering for the last question, wouldn't the probability be 0?

My working:
START -> M1 -> M2 -> M3 -> M4 -> M5 -> M6 -> M7 -> END
Probability Vector [A, T, C, G]
M1 : [0.6, 0.4, 0, 0]
M2 : [0.6, 0.4, 0, 0]
M3 : [0.4, 0.4, 0.2, 0]
M4 : [0.6, 0.4, 0, 0]
M5 : [0.4, 0.6, 0, 0]
M6 : [1, 0, 0, 0]
M7 : [0, 0, 0, 1]
To generate CTAAAAG = 0 * 0.4 * 0.4 * 0.6 * 0.4 * 1 * 1 = 0

[Comments from other post migrated as answer]
The parenthetical statement in the first bullet point says without corrections. In which case, if an unseen nucleotide has a probability of zero, the probability of the sequence is zero. If the unseen nucleotides had the probability of NaN, then the sequence has a probability of NaN, which is a more thorough solution.
The correction would be adding a small iota value to the probabilities to be non zero —fudging against the zero outcome you otherwise find, which may be the result of the training being too small. You have 5 training sequences, so the iota would be <1/5. Let's hypothetically assume that if you were to have picked a 6th MSA sequence it would have had the missing value, so you could set the iota at 1/6. The assumption that you may find that element were you to sample twice as many is more statistically sensible, so iota to 1/10.
I used the letter iota, which is common parlance for something small, but iota is a really bad in maths as it looks like a 1 with English handwriting or in print: you can call the fudge value anything.

source: Stack Exchange

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

sort: votes


text:
NMF clustering not giving results
I'm trying to run the Non-Negative Matrix Factorization(NMF) for my Gene Expression dataset which was originally in matrix form. But it throws errors as given below:
expr <- read_excel("./expr.xlsx")
expr<- as.matrix(expr)

#removing the first column containing gene symbols before clustering
expr_rm <- expr[,-1]
class(expr_rm)
[1] "matrix" "array"

So converted these expr_rm to numeric:
expr_rm <- as.numeric(expr_rm)

#removing negative values for running NMF
eps <- .Machine$double.eps
expr_rm[expr_rm<=0] = eps

estimate <- NMF::nmf(expr_rm, 3, nrun = 10, method = "lee", seed = 123)

Error: NMF::nmf - 10/10 fit(s) threw an error.
Error(s) thrown:
run #1: 'x' must be numeric
Timing stopped at: 1.95 0.14 4.4*

cluster <- cbind(names(NMF::predict(estimate)), NMF::predict(estimate))

Error in h(simpleError(msg, call)) :
error in evaluating the argument 'object' in selecting a method for function 'predict': object 'estimate' not found

result <- list(estimate, cluster)

Error: object 'estimate' not found

For clarity, my expr_rm has row indices as 1, 2, 3, 4,.. which indicates different genes and column headings as samples.
Please suggest a solution for this soon!
Additional information (requested from comments):
>summary(expr_rm) - before conversion
 Sample_C44_0001.genes.results Sample_C44_0002.genes.results Sample_C44_0003.genes.results Sample_C44_0005.genes.results
 Length:145   Length:145 Length:145  Length:145  Class :character

 > summary(expr_rm) - after numeric type conversion
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  0.000   1.808   4.264   4.220   6.333  16.246

If in expr_rm is data.matrix() also, it gives the error like:
> median_centered_data<- expr_rm - med
Error in expr_rm - med : non-numeric argument to binary operator

Answer from @maximilian-press, converted from comment:
Sometimes casting to as.numeric() will convert from two-dimensional matrix to one-dimensional vector. Based on the output you put into comments (which would be helpful in the question), I wonder if this has occurred and caused an issue. For data.matrix, I suggest reading docs for that function- it should have returned only numeric, such that your med variable is likely not numeric, or there is a problem with your data. Possibly the tidyverse tibble created by read_excel might be the issue. I don't know enough tidyverse to say what exactly is going on.

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


text:
Given two ADN sequence, can there be multiple optimal alignement? And how to optimally find them?
Given two sequences, I know we can determine optimal alignment of sequence with the Needleman–Wunsch algorithm. But is it possible that two alignments have multiple optimal alignment (I guess it would mean that multiple alignement have the same score)?
And if possible, how would you find them via an algorithm? Just looking for a general idea on how to go about it, I can't see how you could do it beside by checking for all possible combinations.

As Needleman-Wunsch is a Dynamic Programming algorithm you do need to complete the entire matrix to find the optimal solution and it is possible to have multiple best alignments.

The algorithm assigns a score to every possible alignment, and the purpose of the algorithm is to find all possible alignments having the highest score.

https://en.wikipedia.org/wiki/Needleman%E2%80%93Wunsch_algorithm


text:
Error while parsing gene bank file using Biopython
I am trying to extract the protein sequence of specific genes from gene bank like format file obtained from antismash part of which looks like this
LOCUS       scaffold_10            47160 bp    DNA     linear   UNK 01-JAN-1980
DEFINITION  scaffold_10.
ACCESSION   scaffold_10
VERSION     scaffold_10
KEYWORDS    .
SOURCE
  ORGANISM
            .
COMMENT     ##antiSMASH-Data-START##
            Version      :: 6.0.1-25e43a7(changed)
            Run date     :: 2021-11-05 17:39:45
            NOTE: This is a single cluster extracted from a larger record!
            Orig. start  :: 2049841
            Orig. end    :: 2097001
            ##antiSMASH-Data-END##
FEATURES             Location/Qualifiers
     protocluster    1..47160
                     /aStool="rule-based-clusters"
                     /contig_edge="False"
                     /core_location="join{[2069841:2070128](+),
                     [2070175:2070449](+), [2070507:2071973](+),
                     [2072023:2074178](+), [2074232:2074784](+),
                     [2074838:2077001](+)}"
                     /cutoff="20000"
                     /detection_rule="cds(PKS_AT and (PKS_KS or ene_KS or mod_KS
                     or hyb_KS or itr_KS or tra_KS))"
                     /neighbourhood="20000"
                     /product="T1PKS"
                     /protocluster_number="1"
                     /tool="antismash"
     proto_core      join(20001..20287,20335..20608,20667..22132,22183..24337,
                     24392..24943,24998..27160)
                     /aStool="rule-based-clusters"
                     /tool="antismash"
                     /cutoff="20000"
                     /detection_rule="cds(PKS_AT and (PKS_KS or ene_KS or mod_KS
                     or hyb_KS or itr_KS or tra_KS))"
                     /neighbourhood="20000"
                     /product="T1PKS"
                     /protocluster_number="1"
     cand_cluster    1..47160
                     /SMILES="CC(=O)C(=O)O"
                     /candidate_cluster_number="1"
                     /contig_edge="False"
                     /detection_rules="cds(PKS_AT and (PKS_KS or ene_KS or
                     mod_KS or hyb_KS or itr_KS or tra_KS))"
                     /kind="single"
                     /product="T1PKS"
                     /protoclusters="1"
                     /tool="antismash"
     region          1..47160
                     /candidate_cluster_numbers="1"
                     /contig_edge="False"
                     /product="T1PKS"
                     /region_number="1"
                     /rules="cds(PKS_AT and (PKS_KS or ene_KS or mod_KS or
                     hyb_KS or itr_KS or tra_KS))"
                     /tool="antismash"
     CDS             1899..2720
                     /ID="CDS_13110"
                     /NRPS_PKS="Domain: PKS_KR (26-149). E-value: 7.6e-08.
                     Score: 24.2. Matches aSDomain:
                     nrpspksdomains_jgi.p_FusspF11_1_766668_PKS_KR.1"
                     /NRPS_PKS="type: other"
                     /gene="jgi.p_FusspF11_1_766668"
                     /gene_functions="biosynthetic-additional (smcogs)
                     SMCOG1001:short-chain dehydrogenase/reductase SDR (Score:
                     161.6; E-value: 3.9e-49)"
                     /gene_kind="biosynthetic-additional"
                     /phase="0"
                     /source="prediction"
                     /transl_table=1
                     /translation="MTPPPKSALESDPAVGAQRRFSVTGNAVVTGGAGVLGLHACDALL
                     EHGLEGLMILDVNPAQSQGQITSLQNKFPRAKIMALKVDVTDENVVNAAMEETARVLGS
                     IDTLICFVGVVGCVETLDMPVPQWRKILDINTTGSFICAQAAARQMVKRGRGGSIVFVA
                     SISAHRVNYPQPQAAYNVSKSALLMLKSCLAAEWARYGIRTNSISPGYMDTILNEGDGI
                     AEHRKIWAEHNPSGRMGAPSELTGTVVLLASSAGSYINGADIVVDGGGIVL"
     aSDomain        1977..2345
                     /aSDomain="PKS_KR"
                     /aSTool="nrps_pks_domains"
                     /database="nrpspksdomains.hmm"
                     /detection="hmmscan"
                     /domain_id="nrpspksdomains_jgi.p_FusspF11_1_766668_PKS_KR.1
                     "
                     /evalue="7.60E-08"
                     /label="jgi.p_FusspF11_1_766668_PKS_KR.1"
                     /locus_tag="jgi.p_FusspF11_1_766668"
                     /protein_end="149"
                     /protein_start="26"
                     /score="24.2"
                     /specificity="KR activity: inactive"
                     /specificity="KR stereochemistry: C1"
                     /tool="antismash"
                     /translation="AVVTGGAGVLGLHACDALLEHGLEGLMILDVNPAQSQGQITSLQN
                     KFPRAKIMALKVDVTDENVVNAAMEETARVLGSIDTLICFVGVVGCVETLDMPVPQWRK
                     ILDINTTGSFICAQAAARQ"
     aSModule        1977..2345
                     /domains="nrpspksdomains_jgi.p_FusspF11_1_766668_PKS_KR.1"
                     /incomplete
                     /locus_tags="jgi.p_FusspF11_1_766668"
                     /tool="antismash"
                     /type="pks"
     CDS_motif       1980..2042
                     /aSTool="nrps_pks_domains"
                     /database="abmotifs"
                     /detection="hmmscan"
                     /domain_id="nrpspksmotif_jgi.p_FusspF11_1_766668_0001"
                     /evalue="1.90E-05"
                     /label="PKSI-KR_m1"
                     /locus_tag="jgi.p_FusspF11_1_766668"
                     /protein_end="48"
                     /protein_start="27"
                     /score="17.3"
                     /tool="antismash"
                     /translation="VVTGGAGVLGLHACDALLEHG"
     CDS             join(3252..3263,3321..3443,3493..4428)
                     /ID="CDS_13111"
                     /NRPS_PKS="Domain: PKS_ER (28-277). E-value: 6.4e-13.
                     Score: 40.2. Matches aSDomain:
                     nrpspksdomains_jgi.p_FusspF11_1_794272_PKS_ER.1"
                     /NRPS_PKS="type: other"
                     /gene="jgi.p_FusspF11_1_794272"
                     /gene_functions="biosynthetic-additional (smcogs)
                     SMCOG1040:alcohol dehydrogenase (Score: 275.7; E-value:
                     7.6e-84)"
                     /gene_kind="biosynthetic-additional"
                     /phase="0"
                     /source="prediction"
                     /transl_table=1
                     /translation="MSQTNLSCVLYGPGKARFENRPVPSLKDAHDVIIRISYVGVCGSD
                     VHFWTDGGFARKVSEDQPLVMGHEASGIVRSIGPDVTLLKPGDRVAIEPGFSCRRCKQC
                     KDGRYNLCPKMKFAADPPLTQGTLSRFFSIPEDFAYKIPDSLSLEEAVLVEPLAVAVHG
                     IRLAGLEVGQRVLVQGSGTIGLLTAAVAKAYGAKQVYITDVNLDKIKFAKKYLECSAFI
                     PDLGSTPEENAARFKTETGLDDGVDAVIECTGVEASTQTGLLALSAGGVLVQVGLGKPV
                     QAIPIHAMSEKEIVLKTSFRYGPGDYEIALELLESGKVSVRPLISSITPFEKATEAWEK
                     TRKGKGIKNLIRGVQD"

Now from this I want to extract the protein sequence of gens having PKS domain (PKS_ER/PKS_AT/PKS_KS ...). So far I have tried the following biopython code but it is giving me error
from Bio import SeqIO
from Bio.SeqFeature import SeqFeature

gbk = "myfile_predicted_SMclusters.gbk"
fa = "myfile_PKS_protein.fa"
input_handle = open(gbk, "r")
output_handle = open(fa, "w")

for record in SeqIO.parse(input_handle, "genbank"):
    features = [feature for feature in record.features if feature.type == "CDS"]
    for feature in features:
        if feature.qualifiers["NRPS_PKS"][1] == "PKS":
        assert len(seq_feature.qualifiers['translation'])==1
        output_handle.write(">%s %s\n%s\n" % (
            record.description,
            feature.location,
            seq_feature.qualifiers['translation'][0]))

output_handle.close()
input_handle.close()

The error is as follows:
if feature.qualifiers["NRPS_PKS"][1] == "PKS":
KeyError: 'NRPS_PKS'

The link to the sample inputfile is below:
sample my_file.gbk
I'm using python 3.8.8
Please help!!!

I see 2 problems.
First, you need to check if the qualifier your are searching for exists in the given CDS feature.
Second, the qualifier holds more then just the string PKS for example Domain: PKS_AT (537-846). E-value: 1.5e-98.
You could modify your code like this:
feature_type = "CDS"
feature_qualifier = "NRPS_PKS"
target = "PKS"

def has_pks_value(feature, qualifier, part):
    return any(part in target for target in feature.qualifiers[qualifier])

for record in SeqIO.parse(input_handle, "genbank"):
    for feature in record.features:
        if feature.type == feature_type and feature_qualifier in feature.qualifiers.keys() and has_pks_value(feature, feature_qualifier, target) and len(feature.qualifiers['translation'])==1:
            output_handle.write(f">{record.description} {feature.location.start} {feature.location.end} {feature.location.strand} \n {SeqFeature(feature.location).extract(record.seq)}")

I've created some variables which for your string. This might prevent you from making some typos.
Also instead of a list comprehension, I just run over all features and checked with a single if statement if the asked qualifier exists and contains the given value.
I also changed the string formatting to f-string. But I think you should create the Fasta Files using the SeqIO module of Biopython.
As result I got 37 sequences

Update - I added a missing method to check for the correct value of the qualifier

Thank you for your suggestions, But I god the following error now.  output_handle.write(f">{record.description} {feature.location.start} {feature.location.end} {feature.location.strand} \n {SeqFeature(feature.location).extract(record.seq)}")                                                                                ^ SyntaxError: invalid syntax
could you please Help, I recently started using python (last 3 days) just for this purpose only trying to make sense of resources available online.

Could you try run it without the f-string part? Just use your original version for the file writing

It has another error output_handle.write(">%s %s\n%s\n" % (record.description, feature.location, seq_feature.qualifiers['translation'][0])) NameError: name 'seq_feature' is not defined

The variable seq_feature doesn't exist anymore. Use feature instead


text:
How to split genotypes into two groups from a SNP txt file contains -1, 0 and 1 in R?
I have a SNP file(it is 96 obs. of 1178 variables with -1, 1 and 0 and in some places NA) and I need to calculate Jaccard genetic distance from them in R. Is there any way to convert these 3 values to 0 and 1? I tried to change from -1 to 2 before for a different distance, but the Jaccard distance is binary and you cannot assign 2. Is there a way to re-code Jaccard to take into account -1 and give a good distance matrix? (Jaccard must remain binary and contain only 0 and 1). I will need it later to divide the subpopulation into 2 groups.
Here is an example input data file:
https://pbgworks.org/sites/pbgworks.org/files/snp.txt

THis should be pretty simple to achieve, but could you clarify exactly your rules - i.e. which number gets converted to which?


text:
Is there any function in Biopython to convert a DNA sequence from ambiguous to unambiguous?
I have a project to write in Python that requires to write a function that given  degenerate DNA sequence (for example: KKGTACACCAG) sequence and a molecular weight interval, returns a list of all unambiguous sequences represented by the given sequence (encoding based on IUPAC alphabet: http://www.bioinformatics.org/sms/iupac.html).
Is there any function in Biopython to convert a DNA sequence from degenerate to unambiguous base pairs?
Also, Is it possible to use the molecular weight function from Bio.SeqUtils library for an ambiguous DNA?

from Bio.Seq import Seq
coding_dna = Seq("AAYGANAGYCARAGYAAR")
print(coding_dna.translate())

This gives

NXSQSK

So .translate() will read a degenerate code and provide the unambiguous amino acid, except where it is unable to do so e.g. GAN because this is a split between asp and glu. Thus GAR or GAY is fine, but GAN gives X
You would then use itertools by mapping (map) against the triplet code (import Bio.Data.CodonTable) to obtain all permutations. Doesn't work for X admittedly and stop codons a special case (3 possibilities).
Bio.Data.CodonTable.list_possible_proteins would solve the X problem, but stop codons remain a special case.

I seriously did like the Stackoverflow link in the comments, with the double map funtion ... very slick.

I am missing the SO link the one with :  the double map funtion ... , could you add it to your answer

@pippo1980ONSTRIKE I can't do it now, but I'll look at it early next week. Please check the answer below from acvill

cheers, but R is like ancient greek to me, I am following this question though

M__'s approach is good, but, as explained, some subset of the unambiguous DNA sequences will not be recovered because of the ambiguity in reverse-translation due to the degeneracy of the genetic code. I don't know a python approach, but here's an R function that returns a vector of all the possible unambiguous DNA strings:
make_unambiguous <- function(dna) {
  require(tidyverse)
  require(S4Vectors)
  iupac <- tibble(code = c("A", "C", "G", "T",
                           "R", "Y", "S", "W",
                           "K", "M", "B", "D",
                           "H", "V", "N"),
                  base = c("A", "C", "G", "T",
                           "AG", "CT", "GC", "AT",
                           "GT", "AC", "CGT", "AGT",
                           "ACT", "ACG", "ACGT"))

  tibble(dna) |>
    separate_rows(dna, sep = '(?<=.)(?=.)') |>
    left_join(iupac, by = c("dna" = "code")) |>
    pull(base) |>
    str_split("") |>
    expand.grid(stringsAsFactors = FALSE) |>
    unite(col = sequence, sep = "") |>
    as_vector()
}

make_unambiguous(dna = "AAYGANAGYCARAGYAAR")


text:
Arrange ggplot Figure for scRNA-seq data
I have generated a ggplot for 8 single-cell libraries, with the purpose of visualizing the tSNE facet plot by sample, colored by cell type -- with percentages. The best I could get to is this -

however, it looks too crowded, and I also want the cell types to be labeled for each cluster, which makes it look even worse.
How can I improve the data visualization here?
my code:
p6 <- ggplot(tsnePlot, aes(x=tSNE_1, y=tSNE_2, color=RefinedCellTypes)) +
  scale_color_manual(values=refinedCellTypeColors) +
  geom_point(data = tsnePlotAll, mapping = aes(x=tSNE_1, y=tSNE_2), size = 1, alpha = 1, color = "grey75") +
  geom_point(size = 1, alpha = 1) +
  geom_text(data = tsneClusterLabelsRefined, mapping = aes(x = tSNE_1, y = tSNE_2, label = CellType),color="black",size=0.75,hjust=0) +
  geom_label(data = tsneCellTypeFractions, mapping = aes(x = tSNE_1, y = tSNE_2, label = scales::percent(Freq,accuracy = 0.1)),color="black",fill = "white", alpha = 0.5, label.size = 0) +
  facet_wrap(~Sample) +
  theme_bw(base_size = 12) +
  theme(legend.position="left",
        axis.text=element_text(size=12, color = "black"),
        axis.title=element_text(size=10))

Is there any way I can control the number of cols/rows for the figure?

On top of M__'s suggestion to reduce the point alpha, here are a few other things you can implement to make your figures cleaner:
Remove grid lines from the panel backgrounds
theme(panel.grid.major = element_blank(),
      panel.grid.minor = element_blank())

I think a figure legend like the one you have now is a better approach than labeling every cell type in every facet. However, if you still want to go that route, I suggest implementing repulsive labels with the ggrepel package. geom_label_repel() can also be applied to your percentages so they are not overlapping.

Agreed, BTW I like the user icon

Phage are phantastic!

Except for bacteria :-/

Try
p6 + geom_point(alpha = 1/100)

And variations of 1/10 to 1/1000, this will give a better idea of point density.


text:
DESeq2 throwing error while normalizing raw microarray expression data due to presence of negative values
This question was also asked on Biostars
I am trying to download and analyze a miRNA expression dataset from NCBI GEO (GSE25631). I specifically want non-normalized data to perform normalization and my own set of other analyses later on. Accordingly, I downloaded the Series Matrix File (to match the Sample IDs with their phenotype) and the non-normalized data GSE25631_non-normalized_data.txt.gz programmatically using the GEOquery package in R. To normalize the data I'm using the DESeq2 package. When I run the DESeqDataSetFromMatrix function I get this error:

Error in DESeqDataSet(se, design = design, ignoreRank) :
some values in assay are negative

I went back to the non-normalized data and observed that there are a lot of negative values. I need suggestions on how to deal with this issue. Should I replace all such values with NA? I also observed that every probe has a detection p-value associated with it. Should I also filter out the expression values which have a p-value higher than 0.05?
> head(raw_counts)
        gbm-0001 Detection Pval gbm-0011 Detection Pval gbm-0013 Detection Pval gbm-0024 Detection Pval gbm-0058 Detection Pval
        ILMN_3167304   27.875      0.3071961 -24.9375      0.6689599    14.25       0.350027   4.3125      0.4723163     8.25      0.3966338
        ILMN_3168038  211.875       6.42E-05   1.0625      0.4925717    99.25    0.003645935  60.3125      0.1657064    72.25     0.01086376
        ILMN_3167890 13839.88       3.68E-38 23507.06       3.68E-38 22345.25       3.68E-38 20482.31       3.68E-38 21419.25       3.68E-38
        ILMN_3167526    5.875      0.4577177 105.0625     0.03279042     8.25      0.4117529 -29.6875       0.683706    11.25      0.3604082
        ILMN_3167209   82.875     0.06708142 182.0625      0.0007096    60.25     0.05167192   9.3125      0.4403947    72.25     0.01086376

Here's the R script which I used for all the mentioned steps.
library(GEOquery)
library(DESeq2)
library(tidyverse)

data <- getGEO(GEO = "GSE25631",
           destdir = "E:\\GSE25631",
           GSEMatrix = TRUE,
           AnnotGPL = FALSE,
           getGPL = FALSE,
           parseCharacteristics = TRUE)

clindata <- data$GSE25631_series_matrix.txt.gz@phenoData@data

url = "https://www.ncbi.nlm.nih.gov/geo/download/?
acc=GSE25631&format=file&file=GSE25631%5Fnon%2Dnormalized%5Fdata%2Etxt%2Egz"
download.file(url, "E:\\GSE25631\\GSE25631_raw_counts.gz")
raw_counts <-read.delim("E:\\GSE25631\\GSE25631_raw_counts.gz",
                     stringsAsFactors = FALSE, sep = "\t")
colnames(raw_counts) <- raw_counts[5,]
raw_counts <- raw_counts[-(1:5),]
raw_counts <- raw_counts %>% remove_rownames %>% column_to_rownames(var = "ID_REF")
raw_counts <- raw_counts[ , c(TRUE, FALSE)]

colnames(raw_counts) <- clindata$title
rownames(clindata) <- clindata$title
colnames(clindata)[colnames(clindata) == "disease state:ch1"] <- "Condition"
clindata$Condition[clindata$Condition == "primary glioblastoma multiform (brain tumor)"] <- "GBM"
clindata$Condition[clindata$Condition == "normal"] <- "Normal"
clindata$Condition <- as.factor(clindata$Condition)

dds <- DESeqDataSetFromMatrix(countData = raw_counts,
                          colData = clindata,
                          design = ~ Condition)

Answered at biostars: https://www.biostars.org/p/9561238/#9561240

DESeq2 is meant for digital counts (RNA-seq and similar assays). Arrays are not supported, nor is the implemented method meaningful for array intensity values.
Guides on how to handle these data are probably in here:

Thanks for linking to the other answer, but when you do so, please also reproduce the answer here so that it will not depend on the external link not breaking.


text:
Best R program to use for RNA seq analysis
We have an RNA seq taken from a the same tissue in two different developmental stages. We have one gene, with the number of reads from stage A and stage B, and we want to see which genes are differentially expressed.
After performing RNA-seq, I am left with a CSV file containing a list of genes and the amount of reads in two different cell types. We have one list of genes, with the number of reads in Stage A in one column and the number of reads in Stage B in another column. We are looking to see what types of genes are differentially expressed in the two stages.
I have used DESeq2, but I'm not sure if this is the best/correct to analyze what I am looking for.
What is the best package to analyze the data for significant differential gene expression?
Thanks

DESeq2 works for that.  Unless you literally have only two samples.
There is no statistically sound way of analyzing DE genes with literally only two samples.  This is a bad experimental design, suitable for only the roughest of exploratory purposes.  Only the most obvious differences (that you probably knew already) are likely to be borne out in a legitimate experiment.


text:
Displaying BIC and AICs from MSA nucleotide multisequence fasta file
I am initiating a project reproducing MEGA style trees using R to find the best DNA protein models under AIC/BIC. I have read the R AICs and BIC documentation but I'm struggling to implement this.
Could I implement AICs and BICs without having to specify the number of sequences in the fasta file?
My initial code that works the same as MEGA for aligning and building a phylogenetic tree where I want to view AICs and BIC (as per reference picture):
library(seqinr)
library(adegenet)
library(ape)
library(ggtree)
library(DECIPHER)
library(Biostrings)
library(viridis)
library(ggplot2)
library(msa)

mySequenceFile <- system.file ("sequences", "allSequences.fasta", package ="msa")
mySequence <- readDNAStringSet(mySequenceFile)

mySequence

myAlignment <- msa(mySequence, "Muscle")

myAlignment

alignment <- msaConvert(myAlignment, "seqinr::alignment")
distMatrix <- dist.alignment(alignnment, "similarity")
clustering <- hclust(distMatrix)

plot(clustering)

dendrogram <- as.dendrogram(clustering)
phylotree <- as.phylo(clustering)

plot(phylotree, type="radial")

Please ignore some libraries as I previously ran into some "-" deletion issues when I was writing the code

Having looked carefully at the question, I understand the 'gap in the market' but all this code will do is generate a provisional tree. I would suggest looking at protein models first because you need models to perform an AIC/BIC. So far you haven't generated used or described a protein model and without that there's nothing to compare via AIC.
To answer the direct question you should seek to avoid removing the sample size from the AIC and base the work on the equation below. If you have large samples sizes n can be removed (which is what you are wanting).
The work you should read is jModeltest which is here, https://github.com/ddarriba/jmodeltest2.
This originated from Modeltest by David Pasado.
Their code on GitHub will tell you exactly how to implement AIC and BIC. I recommend ditching R in this instance. If done via subprocessing (below) its not hard.
Their work concerns nucleotides and never really hit the mainstream, but was used by some.
The key criticism is like your approach that the tree generated was based on neighbor-joining. The tree will influence the matrix so there is a danger that the ML or Bayesian phylogeny resulting from the model will be biased by neighbour-joining. Secondly, it was not computationally expensive to do this entirely by ML or Bayesian approaches because it's a 4*4 matrix.
Protein phylogeny is very different, obviously the matrix is bigger 20*20 the other issue is that it's a 'static matrix'. There are approaches where they "combine" matrices but I don't know the algorithm. Thus we know exactly how the protein matrices are generated, "combined" needs a clearer explanation and how this is obtained could be a good application for AIC. Thus could provide useful clarity. There is in addition.the potential to generate new matrices and assess them, but the concept of introducing meaningful dynamism into the matrix does have potential and this is absent (except for cryptic combined options). AIC would then assess whether this is helpful. It might not be hugely helpful in traditional static matrices, but would compare them against further hypotheses of mutation and innovative approaches. If used deployed correctly could be cool because it allows you to experiment and ask the question what precisely does this protein matrix mean for my data?
More generally stretching the utility of AIC type metrics is very interesting in 'bleeding edge' machine learning (admittedly its a bit futuristic).

I suspect (could be wrong) jModelTest is coded in Java, an alternative would be to subprocess via Python using relevant compiled code. This was how the original Modeltest was written, essentially driven through the leading maximum likelihood package of the time known as PAUP. Python does implement AIC/BIC via the statsmodel library, scikit-learn (machine learning) or use the equation directly.

k= number of independent variables to build model
L= maximum likelihood estimate of model
.. for small sample sizes (n)
R is not known for subprocessing and I am not aware of any ML/Bayes libraries for phylogeny.
In summary, this is a good project and it is worth investing the time because it has a clear, meaningful output and potential to extend into more innovative experimental approaches. I for one would use established  code (i.e. written by someone else) that implemented AIC/BIC for proteins matrices and would be fairly straight forward (for you) to write. How you maximise its use in deploying it to help understand protein matrices needs creative thinking (remember its parameters verse likelihood).
Good luck.

Maybe instead of using ML I can try to play around DL as well with some good datasets. Thanks for the the guidance, some things make more sense now and I will try to do everything using python instead.

Cool, DL ... deep learning? Scikit-learn does machine learning, not quite the same. ML/ AIC looks a good staging post (a reviewer will ask for BIC I'm afraid). Please do remember to upvote.

Yes I meant deep learning. I played around some python scripts last night and R as well. Will be creating a GitHub repository later on this week with my progress. Unfortunately I can't upvote yet due to my reputation but I choose it as the answer. Thanks again.

Deep learning isn't trivial and you'll likely encounter skepticism, ML sci-kit lean is widely accepted and much easier.


text:
getting PDB file of interested sequence
I want a PDB file of DNA sequence I am interested in (10 base pairs) for docking. I will appreciate if any of you can please clarify some of my doubts.
Do you think, I can generate PDB file by using this link http://www.scfbio-iitd.res.in/software/drugdesign/bdna.jsp
I did generate the PDB file using the above site. However, when I am using that PDB file for docking it do not give accurate results. It just fits to the sequence and it lacks the information the PDB file downloaded from protein data bank has.
How can I get the pdb file of the DNA sequence I want from protein data bank, I did search it but did not find the one I want. Is it a way to extract the PDB file?
What are other methods I can use to generate PDB file of B-DNA sequence?
Thank you for your time.

Does this answer your question? Generating PDB file for custom DNA sequence

Partially not fully. I already did what was suggested by David.


text:
Downloading exon-wise gene sequence using commandline
I want to download exon-wise gene sequence of around 200 genes from NCBI. Is there anyway to do this from terminal? I have the transcript id for each sequence and I am using efetch command to download their CDS.
eg: efetch -db nuccore -format fasta_cds_na -id NM_000125.4
However, I need their exon-wise sequence as well (i.e. I need the sequence of each exon separately). It will be better if the exon starts with start codon and ends with stop codon (the untranslated regions of those exons completely removed) but I can modify the sequences accordingly once I have all the sequences (UTR included).
Is there any way in which we can give the NCBI transcript id as an input and get exon-wise sequence as an output?


text:
calculate mismatch frequency/rate from a BAM file
I am working with PAR-CLIP data where experimental procedure induces T to C mismatches. The other types of mismatches we can consider coming from artifact such as PCR errors. It could also be an existing SNP. For simplicity, let's consider T to C mismatches are experimentally induced and rest of the mismatches are background mismatch rate.
From Illumina NGS seq after performing alignment, for each of he read groups (reads collapsed in a 30-70 base region with minimum 10 read coverage). Up to this point its already done.
I want to compare T to C mismatches vs background mismatches. To do so I need to calculate the rate of signal mismatches (T to C) and background mismatches.
Using samtools mpileup I can get total number of all kind of mismatches for a read group, however, takes some time.

How do I calculate background mismatch rate from this?
How do I calculate total mismatches in a bam?

I was thinking about something like
for a read group:
mismatch_rate_group=total_mismatches_group/(total_read_group*total_base_group)
for a BAM file:
mismatch_rate_global=total_mismatches_bam/(total_reads_aligned_bam*genome_size_base)

What important factor I might missing in my approach?

I  am not using any variant caller (e.g. bcftools). Initially I tried with samtools mpileup and then used a shell script to figure out the mismatches. But I am not sure is this is the right thing to do. Or even there are any easier methods. If one read has a mutation, I consider that a mismatch and I presume it to be artifact and/or PCR errors.

Please [edit] your question and explain this. You really should use a dedicated variant caller. NGS is error prone, and you cannot just assume that >1 read means the variant is real. In any case, writing your own thing to do this would be far more work and would almost certainly result in worse data than if you use a dedicated variant caller which has been designed for the job. If [edit] your question and give some more details (what sequencing technology, what species, what coverage etc) we will be able to recommend some tools.

@terdon , thanks for your suggestion. I tried to add some more details in my question. I hope this is better now.
Additionally, NGS is error prone : I am less interested in the variant, rather I want to know the rate of such errors so that I can differentiate my signal from all of such errors.

I am not using any variant caller

Use a variant caller. bcftools, at least, will output variant frequency per variant:
bcftools mpileup -ugf ref.fa sample.bam | bcftools call -mv > output.vcf

As an alternative approach, I wrote some code to parse the mpileup output and extract counts for different alleles; this output may be easier for you to parse and calculate your desired frequencies:
# readstomper.pl -- generate base-packed proportion statistics for
#   output from 'samtools mpileup'
#
# example use:
# samtools mpileup -B -Q 0 -f circ-Nb-ec3-mtDNA.fasta LAST_OL_132394_vs_mtDNA.bam | \
#   /bioinf/scripts/readstomper.pl > LAST_OLmpileupProp_132394_vs_mtDNA.txt.gz
#
# example output:
# Assembly,Position,Coverage,ref,cR,pR,A,C,G,T,d,i,InsMode
# Nb_mtDNA,11456,86,T,79.00,0.92,0.0000,0.0465,0.0000,0.0000,0.03,0.01,A;1.00
# Nb_mtDNA,11457,86,A,81.00,0.94,0.0000,0.0000,0.0233,0.0000,0.03,0.02,TA;0.50
# Nb_mtDNA,11458,86,G,79.00,0.92,0.0349,0.0233,0.0000,0.0116,0.01,0.00
# Nb_mtDNA,11459,86,A,78.00,0.92,0.0000,0.0000,0.0235,0.0471,0.01,0.00
# Nb_mtDNA,11460,86,T,5.00,0.06,0.0233,0.8140,0.0233,0.0000,0.08,0.00
# Nb_mtDNA,11461,86,G,57.00,0.66,0.1512,0.1163,0.0000,0.0465,0.02,0.01,AATAACAACACGTAACCGA;1.00
# Nb_mtDNA,11462,86,G,78.00,0.91,0.0349,0.0349,0.0000,0.0116,0.01,0.00
# Nb_mtDNA,11463,86,G,66.00,0.77,0.0465,0.0698,0.0000,0.0116,0.10,0.01,GCAATA;1.00
# Nb_mtDNA,11464,86,T,74.00,0.86,0.0000,0.0814,0.0233,0.0000,0.03,0.00


text:
In Seurat, how can I export cluster IDs to csv files
I am doing scRNAseq analysis with Seurat.
I clustered the cells using the FindClusters() function.
What I want to do is to export information about which cells belong to which clusters to a CSV file.
In a Seurat object, we can show the cluster IDs by using Idents(・), but I have no idea how to export this to CSV files.
I would be grateful if you could show this by using the PMBC data (https://satijalab.org/seurat/articles/pbmc3k_tutorial.html) as an example.

The clusters information, together with additional information on each cells (%mito, sample name, group, orig.ident,...) are stored inside seuratObject@meta.data.
write.csv(seuratObject@meta.data,"./seurat_metadata.csv") will dump everything to a file.
The method @osmoc suggested, is using built in seurat functions to interact with seuratObject@meta.data.
seuratObject$graph_name_resolution points to seuratObject@meta.data$graph_name_resolution.
I personally prefer to use the metadata table, as I have more control on how I interact with the data.

Thank you very much !

Using the meta.data slot is the best answer. The behavior of various shortcut methods can change between different versions of Seurat and even based on how you interact with your Seurat object. More generally, it is a bad idea to take columns from different sources when constructing a data frame.

I've recently done this with my own Seurat data; here's the process:
# load tidyverse packages
library(tidyverse)
# create cluster ID table
tibble(cellID = colnames(pbmc), clusterID = Idents(pbmc)) %>%
  # write out to a file, with today's date
  write_csv(file = sprintf("cluster_mappings_%s.csv", Sys.Date()))

Thank you very much !!

The clusters are saved in the field object$seurat_clusters where object represents your Seurat object.
You can also get them running the command Idents(object).
NOTE: if you run the command FindClusters on different SNN graphs, or different resolutions, you need to specify the clusters by object$graph_name_resolution.
Every time you run FindClusters, the idents of the cells are updated to the latest results, same for seurat_clusters.
For instance, let's say that you have run the command FindNeighbors on the RNA assay (default value).
Then you run FindClusters with graph = "RNA_snn" and resolution=1.
In that case, the clusters are in object$RNA_snn_1.
Once you have the clusters, you can save them as a CSV:
clusters <-  <here you use one of the described methods>
# data.frame
clust.df <- data.frame(cellID = colnames(object), clusters = clusters)
write.csv(clust.df, <path.of.the.csv.file>)

Thank you very much !


text:
KeyError '/n' in Bio SeqUtils Package ProtParam Script
I am trying to run a function from the ProtParam script in the Bio.SeqUtils package. The function is for calculating the instability index of a protein. The code I am using in my script with the pre-built function is the following:
# entry - user inputted string taking in amino acid letters
entry = e1.get()
prot_data = ProteinAnalysis(entry)
result = prot_data.instability_index()

When I run the code, I get the following key error related to the instability_index function.
File "/Genome_Toolkit/venv/lib/python3.9/site- packages/Bio/SeqUtils/ProtParam.py", line 152, in instability_index
dipeptide_value = index[this][next]
KeyError: '\n'

The instability_index() function in the ProtParam script looks like this:
def instability_index(self):
    index = ProtParamData.DIWV
    score = 0.0

    for i in range(self.length - 1):
        this, next = self.sequence[i: i + 2]
        dipeptide_value = index[this][next]
        score += dipeptide_value

    return (10.0 / self.length) * score

I looked at the script myself, specifically the line the error is referring to but couldn't pinpoint the issue. Any guidance is appreciated. Thanks!
the output of print(entry) and print(repr(entry)) looks like this, which is exactly what the entry was in the input:
MALWMRLLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGFFYTPKTRREAEDLQVGQVELGG
GPGAGSLQPLALEGSLQKRGIVEQCCTSICSLYQLENYCN
'MALWMRLLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGFFYTPKTRREAEDLQVGQVELGG\nGPGAGSLQPLALEGSLQKRGIVEQCCTSICSLYQLENYCN'

When the user is pushing enter on their keyboard, it is likely inserting a newline character (\n) at the end of the user input entry and causing the error you have included in your question. You can strip the \n by using entry.rstrip('\n'):
# entry - user inputted string taking in amino acid letters
entry = e1.get()
entry = entry.rstrip('\n')
prot_data = ProteinAnalysis(entry)
result = prot_data.instability_index()

From the comments where the print(entry) was given as

MALWMRLLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGFFYTPKTRREAEDLQVGQVELGG GPGAGSLQPLALEGSLQKRGIVEQCCTSICSLYQLENYCN

and print(repr(entry)) is given as

'MALWMRLLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGFFYTPKTRREAEDLQVGQVELGG\nGPGAGSLQPLALEGSLQKRGIVEQCCTSICSLYQLENYCN'

the \n character is not at the end of the line so .rstrip('\n') cannot be used (it only works at the end of the string).
Try instead .replace('\n', ''):
# entry - user inputted string taking in amino acid letters
entry = e1.get()
entry = entry.replace('\n', '')
prot_data = ProteinAnalysis(entry)
result = prot_data.instability_index()

which on my system (removing the prot_data and result lines) outputs without the \n character and print(repr(entry)) is the same as print(entry).

@Z_Gonzalez rightttt, it's in the middle of the string, I just assumed as it was user input the \n was at the end, .rstrip() only works from the end of the string, we need to use .replace() instead, will edit the answer shortly


text:
What is the limit of partition incompleteness in Bayesian MCMC?
I am seeking advice on whether or not to incorporate incomplete data into my BEAST analysis.
I have HIV data consisting of two partitions of the pol gene: one that is 1.2 kb in length and present in all isolates, & another that is 2 kb in length but only present in one-third of the isolates. I am wondering if I should include the second partition in my analysis, despite its incomplete nature.
I have found some publications discussing the use of incomplete data in phylogenetic analyses, such as Yuchi Zheng & John J. Wiens (2015), Rasmus Hovmöller et al. (2013), John J. Wiens & Jonathan Tiu (2012). However, I remain uncertain about the best course of action.
If you have experience working with incomplete data in BEAST analyses or any insights on this matter, I would greatly appreciate your advice.
Thank you in advance for your help.
Best regards,

Update
As seen from the answer & comments, there has been a misinterpretation of my partitions' structure & their incompleteness. So hereinbelow is a short description of my data, the subject topic, & a diagram.
There are 1009 isolates of the pol gene of HIV from the Republic of Armenia. They were sequenced from the start of the gene by two means: Sanger & an NGS; which caused two groups of length: 1.2 kb & > 3 kb, respectively.
This way the situation may be demonstrated via the next artificial miniature diagram:

Based on the picture like above, I have thought it rational to divide the data into two staple partitions: the first third of pol & its last part. The first partition has all sequences completed, and the second one has mostly full-length complete sequences for a third of isolates and just NNN…NNNs for the other two-thirds.
Ns are automatically added at the alignment step. Sanger sequences are just short and end from 1200 to 1500. I currently do not have access to the source data. They were provided to me already built and aligned at the reference.
In summary, the question is, "Should only the complete first partition be used in subsequent BEAST analysis or should the incomplete second one be included as well?"

@gringer I have updated the question.

Thanks, that's a great improvement; much clearer now.

The original answer was "criticised" for being obtuse.
The basic very simple answer is:

If the data you present in the diagram (alignment) is left EXACTLY in the format used to make the tree - it honestly does not matter whether you have 1 partition or 2 partitions. It makes no practical difference. The posterior likelihood and stuff will change a bit - basically the result is the same. The reason is given in the very long explanation below. I admit 2-partitions looks cooler and there's reason not to do it, computationally it will be a bit slower.
Partitions I am happy to explain why a partition is used  - but thats different answer or question.

What must not be done with HIV is changing every '-' to '?' or changing every '-' to 'N'. Whether this is 1-partition or 2-partitions. With higher eukaryotes that is fine and that will work very nicely, but HIV no. So Weins et al (2012) can do this (almost certainly have done), but you can't.

I know it seems a very small change replacing every '-' with a '?' or 'N', but what it represents seismic. I have attempted to explain why it is so important below the line, but it's about indels (insertions/deletions). I recognise my explanation may not be perfect.
Confusion The confusing bit is you mention 'NNNNNs' and thats the worry. If NNNNs are changed to '-' thats fine, if they are left as 'NNNNs' that raises the complicated spectre of 'indels' and will get very complicated. It probably can be done but its really complicated.
To make the point clear if I perform a trimal analysis on my alignment
trimal -in infile -out outfile -nogap -fasta

It will remove all sites with '-' in them, i.e. it will remove the entire site, not just the '-' character. I think it's really cool. Anyway, the tree from a trimal alignment will be basically the same as '1-partition' or '2-partitions' with all the '-' retained, even though the alignments look very different. That sounds weird but it is a consequence of how Bayes (i.e. Beast) or ML calculate phylogenetic information, which I tried to explain below. HOWEVER, if there is an alignment with 'NNNNs' instead of '-' that will produce a very different tree.
So for simplicity the answer is that 2-partitions is fine for '-', but not if those '-' are represented as 'NNNNN'. Thus for '-' alignments it looks cool, but makes little difference. 2-partitions can make a HUGE difference with different data particularly with loads and loads of NNNNNs instead of '-'.

Summary HIV is a special case because env in particular - and possibly all the HIV genes - are known to contain extensive indels. Indels in env are so prevalent in HIV they can be approximated via distribution (did Ziheng Yang do that?). So the short answer, Weins et al (2012) can do it because they're dealing with higher eukaryotes and genes with negligible indels (in addition they could remove known indel sites), in HIV you would need independent evidence the missing data site was not going to be an indel (i.e. that particular site has never had a known indel in any HIV data set).
The crux is there is a difference between missing data and absent data. Absent data can be a no-go (below; context depending - some times it makes no difference [immediately below]), it can be detrimental to the phylogenetic signal. Missing data is very interesting.
Basically, if you are going to leave the missing 2kb gene as '-' ... thats perfectly fine. It will not add any phylogenetic information - but it will not hurt the calculation. If 'N' or '?' are used instead to represent a missing 2Kb locus for HIV instead  of '-' thats different.

Weins et al (2012) is a very interesting study.
There are four issues

Sample skews
Absent data which causes blockwise deletions in Bayes and ML
Missing data which might recover phylogenetic information at a given site if correctly annotated.
4.Quality of data

The traditional view is that it is better to remove the site/ or sample, if the quality of data is in doubt or for example an entire gene is missing. The traditional view is it is better to be certain of the data regardless of possible risk of a sample skew. However, sample skews are the Achilles heal of population genetics and a lesser extent to phylogenetics. Thats what Weins et al (2012) is concluding, it's really cool. I like it.
We all know that an analyst can't deliberately cause a sample skew (i.e. you can't omit samples without good reason) - that goes with out saying. Right?
The context of the question is "if a sample has pol and its other gene wasn't sequenced - do you fill in the absent gene with "Ns" or "?" to trap the phylogenetic information in pol and the samples which have data in the gene/locus/partition with limited coverage. It's a really good question. My suspicion for HIV is I'm afraid you probably can't, but that might just be env.
Absent data
Bayesian and ML have no '-' character in their matrix. If an indel is present it will trigger that site to be excluded from carrying any phylogenetic information specifically for Bayes and ML. Thus
                Site1
Sample 1         A
Sample 2         A
Sample 3         A
....
Sample 9999      -
Sample 10000     A

Site 1 will not carry any phylogenetic information for Bayes or ML because the single '-' is not represented in the mutation matrix. The entire site for 9999 nucleotides is excluded. A pairwise analysis would be needed, e.g. distance methods.
Missing data However "Sample 9999, Site 1" could be represented as 'N' instead of a '-' because that means any nucleotide in the matrix. That works! Site 1 is now phylogenetically informative under Bayes and ML. Hooray! However, are we certain Sample 9999, Site 1 is a nucleotide rather than an indel? If indels are rewritten as nucleotides, i.e. an absent nucleotide is converted to a nucleotide of unknown type ('AGCT'), i.e. missing data - non-nucleotide to nucleotide - thats cheating.
This creates issues particularly for HIV env where there are loads of indels. If an indel was replaced with an 'N' - the analysis might work nicely, but there is real risk it is not being true to the data.

Possibly explanation Sample skews in phylogeny are really interesting. What can happen is long internal branches can cause changes to the branching order. However, if those internal branches are shorten by more rigorous taxon sampling then these shortened internal branches can recover the true true. I think it was Jo Felsenstein that discovered this (the guy who invented maximum likelihood as a method for phylogenetics in place of parsimony).
One explanation of Weins et al (2012) is that the 'missing data' samples are shortening the distance between internal branches - by providing a more rigorous sample - and that is resulting in recovery of the correct tree. So it is probably a result inline with Jo Felsensteins original observation.
So if pol never has indels - ever - you can do it. This would also help population dynamics and all shorts of calculations. If there's a risk that you are replacing indels with nucleotides it is very much, much better to avoid it.
Conclusion Everyone in phylogentics knows the issues of sample skews - a possible sample skew is permitted if the data is genuinely never been produced - it's not in existence (obviously). However, it is not permitted if the data exists but the analyst hasn't included it - for whatever reason. Its also preferred if the quality is in question - its better to err on the side of caution.
If you fill in an absent gene with "?" or "N" to obtain phylogenetic information from genes that are present for that sample ... that is permitted providing it doesn't change the possible nature of the underlying data. The rationale for doing this is the risk of sample skews.

So would I adopt a "non-traditional" approach? As a general rule for viruses no and this is a particular special case regarding indels.
If this was a higher eukaryote and like - that critical sample and  that gene - didn't give sequence. This has happened. So would I dump the sample (or in this case the locus) - no I'd caveat the tree and explain and empirically check the assumptions.
If I was in a situation where 2/3 of the data were missing, for a higher eukaryote (DNA with great proof-reading), personally I would do a full scale investigation, i.e. phylogenies with and without the gene/locus/partition in question. That is pretty much was Weins et al (2012) have done.

Thank you very much for your reasoning! My question was mostly dedicated to the whole partitions, not single incomplete sites or indels. My sequences jointly cover only pol, but it also has a bunch of indels. I am intent on not including insertions and exchanging deletions with N.

@M__ I am really thankful for your attention! But please, could you shed light on your final point about indels? Have I understood you correctly that we should exclude all insertion columns and exchange deletions with N (I was intent on doing this initially), shouldn't we? Or we should exclude deletion sites as well, or we cannot analyse HIV fragments susceptible to indels with Bayesian MCMC at all?

@Vovin I've summarised the basic answer above the question. Please check it and if it still not clear I will attempt another explanation.


text:
Is it normal to get truncated identifiers after applying SQL scema of RegaDB?
This question was also asked on Biostars
I am trying to set up RegaDB on my CentOS system, and I have installed PostgreSQL via conda. I followed the RegaDB installation manual to configure the PostgreSQL database. However, when I have applied the postgresSchema.sql schema from the RegaDB software package, I noticed several NOTICE messages about truncated identifiers.
Here is a sample of the output I received:
psql:postgresSchema.sql:25: NOTICE:  identifier "query_definition_parameter_type_query_definition_parameter_type_ii_seq" will be truncated to "query_definition_parameter_type_query_definition_parameter_type"
psql:postgresSchema.sql:27: NOTICE:  identifier "query_definition_run_parameter_query_definition_run_parameter_ii_seq" will be truncated to "query_definition_run_parameter_query_definition_run_parameter_i"
psql:postgresSchema.sql:29: NOTICE:  identifier "resistance_interpretation_template_resistance_interpretation_template_ii_seq" will be truncated to "resistance_interpretation_template_resistance_interpretation_te"

I understand that PostgreSQL has a maximum identifier length of 63 bytes, and identifiers exceeding this limit will be truncated. My question is whether these truncated identifiers will cause any issues with the proper functioning of RegaDB. If so, what would be the best approach to fix this issue?
I would appreciate any guidance on this matter, and if there are any additional precautions or testing measures I should take to ensure RegaDB functions as expected.
Best regards,

Hi @Vovin I would make sure RegaDB is not abandonware. Its github hasn't been updated for 9 years. PostgreSQL will have been revised many times since then and will not be the version used to develop RegaDB. Upvoted anyway.

@M__ I suppose that newer versions of PostgreSQL (as well as Tomcat) actually caused a bunch of my fail attempts to run RegaDB. But the maximum length of identifiers is a substantially base thing in PostgreSQL, as it is set by the constant NAMEDATALEN at the compilation step.


text:
How to Sort and Index a SAM file without converting it to BAM?
Generally, I use samtools to sort and index BAM files. Samtools works also with BGZF compressed SAM but not with plain ones.
I have terabytes of plain SAMs at the moment and I would love to avoid compression/conversion. I have all the data in several external HD and I am not allowed to delete the original SAM files. I would need much more space than I have now to also store the BAM files.
So, I was wondering if there is a way to avoid it. Is there some tool/trick to sort and index uncompressed SAM files, without converting them to BAM or compressing them?

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

In general, avoid storing any SAM files on disk. Instead, use htslib/pysam/etc to directly read BAM files, or pipe SAM to other tools or your own scripts.

Not only will you save disk space by converting to BAM, but BAM files are faster to manipulate than SAM. Source: Dave Tang's SAMTools wiki.
sort supports uncompressed SAM format from a file or stdin, though index requires BGZIP-compressed SAM or BAM. I don't think you can get around this.
If you're dead-set on reading from and writing to uncompressed SAM format, you can use pipes to avoid writing BAM files. Specify the output of view to be an uncompressed bam, then pipe it to your next command, then pipe it to view again to output your sorted SAM file.
samtools view -u in.sam | samtools sort | samtools view -h > sort.sam
samtools view -u sort.sam | samtools index

I suspect you could convince tabix to index a SAM file, though you're correct that using BAM is a MUCH better solution.


text:
Allele effects by parental-origin estimated from genotyping array data
We are dealing with relatively large family data (> 1,000 trios) with genotyping array data.
We would like to see the allele effects on children separately of mothers, fathers, and children, as in a previous study (PMID: 34282336).
Unfortunately, this is difficult for us to replicate, so we need another method.
Is there any software that can achieve this or a similar traditional approach?
Also, can that be applied to grandparents?
In my understanding of the above mentioned study is they:

identified the parental genomes within the child genome;
thereon conducted a GWAS for birth weight

Overall, the aim of the previous study was to estimate the alleles responsible for birth weight or in our situation a comparable phenotype. It is important to note the prior study did not use publicly available algorithms.
Thus, we wish to identify the maternal and paternal alleles within their progeny via public domain algorithms, and therefore do not necessarily need to exactly reproduce the previous study.
We already have imputed data, and are using GCTA for usual GWAS. To achieve our goal, we need to identify the transmitted and non-transmitted alleles in parents.

Hi @PenguinPartyH0, thanks for the extra info. Just to be clear I don't do eukaryotic genomics in general, but I've edited the question with your added information. It would be useful to explain what commercial algorithms the authors used, because it is likely someone here will know the freeware version. That is very likely the most important information required. You can include this info by editing the question. I personally don't think this is a problematic situation. I have never personally known any commercial algorithms that are more powerful than those in the public domain

Has your data been phased as well as imputed, or just imputed?

We have just imputed data and genotyped data.

The first haplotyping step is typically done in two stages:

Imputation of missing genotypes - can be done using the Michigan Imputation Server
Phasing of called genotypes to create semi-contiguous haplotypes - can be done using SHAPEIT

I have many thoughts about the best approach for GWAS (mostly around the idea of population subsampling), which could potentially also apply to what you want to do. Unfortunately, it's been a long time since I did family population genetics, and I've forgotten the things we did to deal with complex population structure. I'm pretty sure we used PLINK for at least some of it, but I notice that PLINK has a disclaimer about close family structure, which suggests it may not work so well for trio data:

This method does not properly adjust for small-scale family structure. As a consequence, it is usually necessary to prune close relations with e.g. --king-cutoff before using --glm for genome-wide association analysis. (Note that biobank data usually comes with a relationship-pruned sample ID list; you can use --keep on that list, instead of performing your own expensive --king-cutoff run.) If this throws out more samples than you'd like, consider using mixed model association software such as SAIGE, BOLT-LMM, GCTA, or FaST-LMM instead; or regenie's whole genome regression.


text:
Using Diffusion Map on normalized - or PCA reduced data?
This question was also asked on GitHub
In the R package destiny, a Diffusion Map is often performed on the first 50 PCs of the data.

Is there a theoretical reason for performing a Diffusion Map on PCs rather than on the data itself ( other than speed)?

Is it reasonable to do the diffusion map on the normalized data itself rather than on Principal Components of the data, if that differentiates the states better ?

Techinically, I am using the  R Bioconductor package destiny. The help page on destiny::DiffusionMap tells "Using e.g. 50 PCs results in more regular looking diffusion maps", but it's not clear compared to what will they be more regularly looking.

The empirical statistical point of view
The technique is  uniform manifold approximation and projection (UMAP) and is comparable to t-SNE. UMAP is considered a slight improvement, personally I think they have similar behaviour. They are dimensionality reduction methods.
The reason UMAP or t-SNE is often not applied directly is it is weaker than PCA are capturing the multi-dimensional variance of the data. Both methods can be applied directly to the data, but will not necessarily deliver stronger clusters than PCA (which is the goal of the technique).
However, the weakness of PCA is that the variance can be spread across 50 principle components, making it hard to recover the representative clusters in certain situations. These situations are commonly found in RNA expression analytics, notably where the first 3 or 4 PCs only comprise ~30% of the variance. PCA should be trapping >80% of the variance within the first 3 or 4 PCs. In expression data (RNAseq) the next 46 PCs might comprise a total of ~30% of the variance and omitting that is mistake, but manually comparing 50 PCs in bivariate plots is not practical. So at best a good unsupervised learning algorithm (UMAP/t-SNE) following PCA in RNA expression via RNAseq might trap 60% of the variance. It's not ideal but it's okay. Thus you can't ignore the 30% of the variance in the bottom 46 PCs and UMAP/t-SNE is recovering that. If 80% of the variance was in the top 4 PCs, then the remaining 20% of the variance in the lower 46 PCs can be ignored and t-SNE or UMAP is not needed.
Resolution The way you can capture the variance from all 50 PCs is to perform a clustering method over the top of PCA, this being t-SNE or UMAP.
Summary,

t-SNE or UMAP without PCA might trap at best 20% of the variance.
The first 4 PCs in a PCA might trap 30% of the variance (it will be a better result than t-SNE or UMAP).
Combining PCA and then UMAP/t-SNE will capture 60% of the variance.

60% of the variance, its just about acceptable. Thus thats why this approach  is done from an empirical point of view.
Practical point of view
PCA and then UMAP (or t-SNE) will give tighter, neater clusters (cleaner, less messy groups), but thats not the reason why you are using it.

The publications in the link,
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC10312943/
and the PLoS Biology paper
https://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.3001532
... both use UMAP across PCA. Its a major approach to cluster resolution.
Diffusion maps are an attempt to extrapolate a variance within Eigenvalues via a random walk, i.e. to defuse within Euclidean space. It is not a primary calculation in my personal opinion because you are sort of introducing more variance into a variance dimensionality reduction calculation. Once you have very tight clusters using PCA/UMAP or PCA/tSNE then you could check to see whether Diffusion map clusters are similar.
This is my perspective only, not to replace PCA with Diffusion map nor to consider a "down-stream" calculation. Its separate calculation used to verify if PCA-UMAP is okay.
If you doing graph theory a Diffusion map could be really important as a null hypothesis and a really important way to obtain a different understanding of the graph (i.e. this is a technical term for a network). Using a diffusion map without any directionality (the papers you cited)*, I struggle to see this as a primary analysis. In the PLOS Biology paper is part of the appendix.
*, graph says A->B->C
The publications treat A,B, C as unconnected. If you can hypotheses how A, B and C are connected then diffusion map could be cool. Thats my opinion from the perspective of graph theory.
Summary I would definitely choose PCA-UMAP over a diffusion map every time. This is what has happened in the two publications cited. The diffusion map is secondary analysis, so by all means do it ... but present PCA-UMAP first, and base the results around that.

From the comments there's some sort of time-series. So if you are performing a time-series - it's not datamining - this would be the justification of using a Diffusion map over PCA. This therefore answers points 1 and 2. However, it would imply you are pooling your data.

The specific technique is a Diffusion Map after PCA (and not UMAP/t-SNE). Does the same reasoning hold for the Diffusion Map dimensionality reduction as well?

Thanks for clarification, I've answered above

I'll explain my specific case in more detail : The reason I personally use a Diffusion Map is not only to see the clusters, but to draw a trajectory between the clusters. In my case, not only the Diffusion Map separates better between the (two) main clusters, but the genes change in a meaningful way, across the trajectory : from state A, to intermediate state, to state B.  However, this is not stable in the sense that it appears only when using the normalized data as a basis for Diffusion Map, not when using the 50 PCs.Hence, the two OP questions focus on Diffusion Map(on normalized vs on PCs)

@Sam maybe this is a separate question? The studies you presented suggested a straight cluster analysis. If you are doing e.g. a time series that will require a lot more thought.

I have intended the studies just to illustrate the point that taking the first 50 PCs is often done (so that I do not get called out for making unfounded statements about what is often done). They were not intended to be the focus of the question, but I see it was not clear.
How would you phrase a separate question? In my mind, the question is well-phrased (but it might be I don't get it).

Are you doing a time-series?

My data has been collected at one time-point, but computationally I perform trajectory analysis using Slingshot, based on that, divide the trajectory into 3 timepoints (start, mid, end), and perform differential expression test between the timepoints using tradeSeq

can you help further?

Thanks @Sam answered above