text:
Analysis of 4665 proteins in String database
Background I have a data set of 4665 proteins from Phytophthora cactorum and want to analyze them in the String database as other database does not have Phytophthora cactorum for GO or KEGG analysis. The String database states I can upload a maximum of 2000 proteins at once.
Thus I subset my data into 3 sets and upload. However, it gives/risks statistically biased analysis.
Question Is there any possibility I can analyze all of my 4665 protein sequences at once?
Looking forward for the help.
Best regards,
Bikal

I'll migrate this to bioinformatics which it is better suited.

What kind of analysis? Are you looking for the proteins each of your 4665  interact with? Something else? Please [edit] your question and explain the kind of analysis you want to make, the questions you will be asking, so we can understand and suggest alternatives.

Do you just want to annotate them with function, or something else?

@terdon I wanted to have GO and KEGG terms and subcellular localization prediction as string website gives when we upload a set of protein sequences. SO my objective is to have GO and KEGG enrichment analysis. I do not need network. No other database has Phytophthora cactorum for GO and KEGG analysis.

I didn't even know STRING has the kind of analysis you are requesting, I don't see it on their webpage. In any case, the standard tool for GO or KEGG analyses, at least back when I was doing this sort of thing, is DAVID.
I have confirmed that it does have Phytophthora cactorum so all you need to do is upload your file there and use their tools.

String is a major protein database. Phytophthora infestans is on KEGG but I am not aware that Phytophthora cactorum is there. I think there could some confusion. For the OP the difference is likely important.

@M__ yes, I know STRING but I knew it as a database of protein protein interactions and wasn't aware that they also offered other tools. But I am not suggesting KEGG, I am suggesting DAVID which was my goto tool for GO enrichment and pathway analysis when I was working in this field. And yes, it does have P. cactorum, I checked before posting.

Okay I don't know DAVID.

Just to confirm the OPs question
going to
https://string-db.org/cgi/organisms and typing "Phytophthora cactorum", then clicking 'random network' results in..
https://string-db.org/cgi/network?taskId=bbh7BVTwU0tH&sessionId=bA71KKVVotqM
I've encountered String for metabolic pathways, so I don't have direct experience with the OPs query.

All databases are all limited by skews in underlying wet-lab work, e.g. there are some large taxonomic groups of bacteria that are complete blanks, i.e. no-one has ever worked on them. This however doesn't concern the OP because they are focussing on a single plant species (I think thats what it is, but could be wrong).
The limit on uploading is understandable - these are not our computers and someone else pays for the OPs analysis.
What the OP suggests is the '3 uploads' comprise a sampling bias and this too is very, very understandable regardless of what analysis is performed.
Example If you've loads of COG 'V' type genes (disease resistance), for example, in one of the uploads, but not evenly mixed this will certainly skew the results.
Solution The solution in this scenario is to randomly sample the genes amongst the 3 uploads. This is easy to do in Python via random and itertools libraries/packages . It would be easy to write a custom script. (separate question).
Normally in a randomisation procedure there needs to be several, e.g. 10, iterations. So essential the OP would need to upload and perform the analysis 30 time comprising 3 groups of 10 replicates.
Example 2 Equally if the OP is honing a particularly metabolic pathway and the genes got split between the three uploads then visibility of the pathway would be partial.
Solution The OP could cluster genes into one of the three uploads by membership to a given metabolic pathway. This is more difficult. This will depend on what the upstream annotation software/algorithm was because often the pathway membership is described in that annotation file. string is a strong database for pathway membership.
Example 3 The skew is taxonomic membership, I don't this this counts here because it is a single species.
Conclusion

Ultimately, the preferred solution is to perform the analysis locally this way the OP will not need to break up their file. With KEGG the databases can be downloaded.
The OP could contact string essentially asking if they could supply an account that lifts the upload limit. string would likely want a fee if they will do this at all. I honestly have no idea.
There are other servers as @terdon recommends, I think BRITE is another one which might have the plant (which is what I think it is) under investigation.
The OP could go through solutions 1 and 2 above. Solution 1 is likely to be 'the solution', but there needs to be a supporting script written (separate question).

APIs I don't know about string but if there was an API this would enable a multiple replicate analysis to be coded without all the manual point and click that it would otherwise involve.

My opinion of Go is that it is sparsely populated, just saying.

Hej Bikal,
I am not sure why you would like to inspect all 4,000+ proteins at once.  I would expect that you have a look on a subset of "interesting" genes, eg from a differential expression experiment.
Yet, since you've asked, yet it is possible to get the data for all proteins at once.
You can download
the network for entire genomes
and analyze with what ever tool you seem fit.
Alternatively, may I recommend you to have a look at the Cytoscape plugin stringApp. It should help you with the download, analysis (including pathway enrichment), and visualizations. Unlike the String network in your web browser, Cytoscape can easily handle large networks.
The tool's authors have a nice series of
video tutorials of how to use it. Hopefully that helps in your work.


text:
Using information criteria weights to create consensus tree
While I am familiar with calculating information criteria for various ml models of a multiple sequence alignment, I am not aware of how to use the derived weights to create a tree that combines multiple models weighted accordingly.
Thus, I have X trees, computed by multiple different models on the same alignment. Each has an associated BIC. BICs can be converted to weights. I wish to average the trees according to the BIC derived weights.
I would prefer to be able to do this with R.

Calculating a consensus is easy, RAxML has a tree consensus module, however calculating it according to weights is a rare procedure. ETE3 (Python) will have a consensus library and I assume in R it will be ape
The only way I can see this being done is to use the relative BIC as a duplicate factor of the number of trees of a given model and then perform a global tree consensus on the resulting treefile. Thus if Bayes Information Criterion top value is twice the size as the lower value the number of 'top value BIC trees' is double that of the lower value. That could be coded up to be cleaner. . Alternatively a multiplication factor could be introduced to ensure the number of duplicated trees within a scale of 100 - to give a cleaner more interpretable 'bootstrap' (but its not a bootstrap) score for each branch.
Here the consensus frequency is a direct function of BIC, essentially you just have an alternative measure to the "non-parametric bootstraps".

I better finish this question. The solution to the problem is to produce separate trees for each model, provide a method of robustness, e.g. bootstrapping and then produce a single tree with all the bootstraps for each model along the branches. That is the formal approach if you're not going to simply select the model with highest probability, highest BIC in this case. This however is not the question.
The problem with the solution above the line ... is the question. I was sort of being polite when I said 'rare procedure'. Weighting a consensus tree is not done in this context, weighting a locus or specific site - alignment position - can be done with good justification within the tree building method. The problems sort of multiple when it is weighted according to exact difference in BIC because the probability of a given branch being robust - which does use consensus is calculated differently. What you are doing is assigning shifts in the BIC/maximum likelihood not to the global topology but to specific incongruent branches. Thats not how it works. The BIC/maximum likelihood - whatever - is for the whole tree. The difference between models is the whole tree - every single branch change - not just a branch here and a branch there because its the score (probability ML) of a singular hypothesis. So a branch change here or there is a different hypothesis - for which you might not have generated a BIC score. The branch here and a branch there is MCMCMC consensus or bootstrap resampling consensus not precise shifts in BIC/maximum likelihood. You would need to dedicate the entire project to this specific method and even then it might be not accepted ... its not as a passing calculation.
There could be biological situations where weighting a consensus tree is doable in the way described here. If a specialist reviewer saw the calculation they would refuse to permit publication because it violates the foundations of tree building theory. What you may not be aware is the long historic that has led to the current acceptance of modern methods and what you are proposing is dismantling that theory.

I think I can make something similar work. Weights are not linearly related to IC or to IC ratios. They are proportional to exp(-delta/2), where delta is the difference between the individual IC and the lowest IC. Thus, the lowest IC automatically has a weight of 1, and each successive IC has a lower weight, therefrom. But I can make that work.

source: Stack Exchange

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

sort: votes


text:
Filter with bcftools
I want to filter a SNP, specifically CHROM1:POS:630128. I tried to use bcftools for this.
I came up with the following variant of a bcftools call:
bcftools filter -e "CHROM=1&POS=63018" /home/oem/Documents/GitHub/Genomics_prac_guide/trio/files/genome_file.vcf.gz

which leads to:
Error: cannot use arithmetic operators to compare strings and numbers

and the call with unquoted parameters leads to
$ bcftools filter -e CHROM=1&POS=63018 /home/oem/Documents/GitHub/Genomics_prac_guide/trio/files/genome_file.vcf.gz

sh: 1: /home/oem/Documents/GitHub/Genomics_prac_guide/trio/file/genome_file.vcf.gz: Permission denied
Failed to read from standard input: unknown file type

Q1: Is bfctools filter the correct tool here to use?
Q2: Regarding this comment the first error might indicate that I have to annotate the CHROM & POS Tags in the header file. Is this correct?

An example of the file:
##fileformat=VCFv4.2
##fileDate=20221104
##source=PLINKv1.90
##contig=<ID=1,length=248928329>
##contig=<ID=2,length=242157857>
##contig=<ID=3,length=198157342>
##contig=<ID=4,length=190106194>
##contig=<ID=5,length=181288140>
##contig=<ID=6,length=170595721>
##contig=<ID=7,length=159331484>
##contig=<ID=8,length=145075470>
##contig=<ID=9,length=138211490>
##contig=<ID=10,length=133659511>
##contig=<ID=11,length=135075227>
##contig=<ID=12,length=133261768>
##contig=<ID=13,length=114341522>
##contig=<ID=14,length=106873816>
##contig=<ID=15,length=101922277>
##contig=<ID=16,length=90172890>
##contig=<ID=17,length=83203771>
##contig=<ID=18,length=80252738>
##contig=<ID=19,length=58586567>
##contig=<ID=20,length=64323573>
##contig=<ID=21,length=46679699>
##contig=<ID=22,length=50776369>
##INFO=<ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on real reference genome">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  FAM001_ID001
1   69869   rs548049170 T   .   .   .   PR  GT  0/0
1   74792   rs13328684  N   .   .   .   PR  GT  ./.
1   630128  rs9283150   A   .   .   .   PR  GT  0/0
1   791532  i713426 N   .   .   .   PR  GT  ./.
1   792461  rs116587930 G   .   .   .   PR  GT  0/0
1   817341  rs3131972   G   .   .   .   PR  GT  0/0
1   818725  rs12184325  C   .   .   .   PR  GT  0/0
1   820888  rs12567639  A   .   .   .   PR  GT  0/0
1   823656  rs114525117 G   .   .   .   PR  GT  0/0
1   841166  rs12124819  G   A   .   .   PR  GT  0/1
1   858952  rs12127425  G   .   .   .   PR  GT  0/0
1   866156  rs79373928  T   .   .   .   PR  GT  0/0
1   880041  rs72888853  N   .   .   .   PR  GT  ./.

Can you give us an example of your file?

Done! See edited post.

First of all, bcftools might be a bit of overkill if you want one specific variant. All you need is GNU grep (the default on Linux):
zgrep -P '^(#|1\t630128\t)' file.vcf.gz

Or, with any grep:
zgrep -E '^(#|1'$'\t''630128'$'\t'')' file.vcf.gz

That will print all lines that start with a # (the headers) or that start with 1, followed by a tab, then 630128 and another tab.
Alternatively, you can use awk:
zcat file.vcf.gz |
    awk -F'\t' '/^#/ || ( $1==1 && $2==630128)'

All that said, yes, you can also do it with bcftools and in most cases that is the right tool for the job so you might as well use it here too to get used to it. However, you need the right command and that's view, not filter:
bcftools view -r 1:630128 /home/oem/Documents/GitHub/Genomics_prac_guide/trio/files/genome_file.vcf.gz

Note that the file needs to be indexed, so if you have not done so already, run
tabix file.vcf.gz

Also, the reason your unquoted command failed as it did is because the & has special meaning to the shell, it means "run this command in the background". Therefore, when you ran bcftools filter -e CHROM=1&POS=63018 file.vcf.gz that actually ran two separate commands:
bcftools filter -e CHROM=1&

And then:
POS=63018 file.vcf.gz

Which means "set the variable POS to the value 63018 and execute the command file.vcf.gz. Since the file isn't executable, you got:
sh: 1: file.vcf.gz: Permission denied

First of all thanks for the quite detailed answer! This gives me a good introduction into bcftools environment. So when I use view to filter SNP's the reason for  -filter is then to apply some standard filter to the file?
Second question: If I want to remove this single variant from the file is the correct way to do this bcftools view -r 1:630128 -e file.vcf.gz ?

I think filter is more to get things like "show me variants where the AB is X and the GT is Y" and things like this, while in this case you are just focusing on a specific region, albeit a region of length 1, so view is the command for this. Reversing the match is a completely different question. I don't know if bcftools can do that, but zcat file.vcf.gz | grep -vP '^1\t630128\t' will get rid of all variants not matching that specific position.

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


text:
Using very closely related strains to increase coverage for short read, de novo assembly
Do you think it's possible to combine short read Illumina libraries (WGS) from multiple closely related eukaryotic microbial strains (e.g. libraries from a re-sequencing study, >99% ITS1 sequence) to boost read coverage during a de novo genome assembly run (e.g. with SPAdes)? My target genome is much larger than anticipated, and I'm wondering if increasing read coverage by an order of magnitude could improve length accuracy, contiguity, etc.
Maybe differences in repetitive regions or transposable element movement prevent such an approach? Or maybe kmer frequencies are affected enough that de brujin graph construction gets messed up.
There seem to be several pitfalls, but maybe benefits could outweigh the costs here. If anyone else has done this or considered it, I'd be interested in your thoughts. Thanks!

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.

Given that the differences between your strains probably involve rearrangements and copy number variations, this approach seems inadvisable if your goal is a highly contiguous genome. You will make the dBG into more of a hairball if such differences exist.
Illumina is quite cheap these days, I would vastly prefer simply sequencing more deeply if your issue in inadequate coverage. Online tools exist for figuring out how to get there.
Unfortunately, you will swiftly plateau in terms of contiguity even for quite simple genomes using Illumina. So you can expect quickly diminishing returns.
Knowing more stats on your assembly would help evaluate this. For example, if you are working with effectively 10X coverage with N50 ~ 500bp, I might be willing to look the other way while you co-assemble such a synthetic metagenome, because you will at least be able to assemble gene length contigs, which are unlikely to vary too dramatically between your strains. However, if the issue is that your assembly has an N50 ~ 10Kbp and you want it to be higher, then you should consider using other approaches to improve your assemblies, such as long reads. In a pinch, you could do something like optical maps or Hi-C, though I believe that is not preferred these days given much cheaper/better long reads. (Full disclosure: I used to work at a Hi-C company and still hold stock.)
The field of metagenomics is based on the assumption that your suggested approach is ok to get a basic idea of e.g. gene complement, and if you are clever you can even deconvolute such mixtures into genomes, under the assumption that different genomes are at most distantly related. But that doesn't mean that it's a good idea to do it in silico by combining isolates.

Thanks for your thoughts. I actually do have PacBio HiFi contigs as well, but hybrid assembly efforts have only gotten me so far (tho much higher N50 than you mention, which is nice). Also not in position to continue sequencing :( I'm sorry if this is too speculative for this forum, I've just never seen anyone publicly entertain this notion. I think i agree that broader rearrangement patterns will likely hurt more than the process might help.

@bishopia I would say that if you have HiFi contigs, go with those! You are not going to improve them by adding a bunch of heterogeneous shotgun Illumina data, as you seem to have inferred.


text:
Large amount of REF prefix differences after lifting
I lifted some unmerged trio files from hg19 to hg38 using liftOverPlink. It looks like that several houndred allels on each chromosome of one file shifted/switched to non-standard alleles. I recorded some of the differences. And will add them to the post. I will also add a sample of the file to this post. My goal is to merge the lifted trio files to one singe file.
My command to merge is
bcftools  merge -l merge.txt --force-samples -Ob -o file.bcf

The REF prefixes differ: C vs N (1,1)
Failed to merge alleles at 1:3940550 in file_1.vcf.gz

First, I tried to remove those alleles in question but I realized that these are just to many to do this by hand.
I used R to edit the file
file <- read_delim(
  gzfile(
    str_glue("{some_path}file.vcf.gz")),
  delim="|",  quote='', escape_backslash = T, col_names = F) %>%
#A list of all the REF prefix differences: number equals position in chromsome 1
  #630128 N vs A, 791532, N vs A, 841166 G vs A, 937572 C vs T, 1110226 N vs A, 1486611 N vs G, 1554290 C vs G,
  #1654663 N vs T, 1671724 C vs T, 1755523 C vs T, 1774697 C vs T, 1916092 C vs T, 1949922 C vs T, 1959595 C vs n
  #1959717 C vs T, 1983048 A vs G, 1991690 A vs G, 1999315 C vs T 2096621 C vs T 2107385 A vs G, 2124676 A vs G
  #2274142 A vs C, 2284952 A vs G, 2293812 C vs T, 2303949 N vs C, 2321888 T vs N, 2351827 C vs T,2364905 N vs G
  #2364905 T vs C, 2420707 T vs C, 2431388 T vs C, 2441501 T vs C, 2461209 C vs G, 2614790 A vs G, 2620461 G vs A
  #2622185 C vs T, 2792599 A vs C, 2861494 C vs T 2864702 G vs N, 2884008 T vs C , 2894837 A vs C, 2925918 T vs C,
  #2944366 T vs C, 3097126 G vs A,, 3107701 C vs A, 3109422  C vs A, 3112423 G vs A, , 3119756 C vs T, 3146311 G vs T
  #3282557 G vs A,, 3294386 A vs G, 3294894 A vs C, 3303337 A vs G, 3320610 T vs C, 3335540 T vs C, 3344952 T vs C
  #3345118 T vs C 3355977 G vs A, 3401871 G vs A, 3443779 G vs T,, 3534441 T vs N, 3591297 A vs G,
  #3595370  A vs G, 3601923 T vs C, 3624792 G vs N, 3689917 T vs C, 3689917 G vs A, 3695328 G vs A, 3696503 G vs A
  #3762340 T vs C, 3775163 A vs C, 3870048 T vs Cl 3924314 G vs A3935847 T vs C, 3940550 G vs A
   filter(!grepl('.+630128.+|.+791532.+|.+841166.+|.+937572.+|.+1110226.+|.+1486611.+|.+1554290.+|.+1654663.+|.+1671724.+|.+1755523.+|.+1774697.+|.+1916092.+|.+1949922.+|.+1959595.+|.+1959717.+|.+1983048.+|.+1991690.+|.+1999315.+|.+2096621.+|.+2107385.+|.+2124676.+|.+2274142.+|.+2284952.+|.+2293812.+|.+2303949.+|.+2321888.+|.+2324778.+|.+2351678.+|.+2351827.+|.+2364905.+|.+2364905.+|.+2420707.+|.+2431388.+|.+2441501.+|.+2461209.+|.+2614790.+|.+2620461.+|.+2622185.+|.+2792599.+|.+2861494.+|.+2864702.+|.+2884008.+|.+2894837.+|.+2925918.+|.+2944366.+|.+3097126.+|.+3107701.+|.+3109422.+|.+3112423.+|.+3119756.+|.+3146311.+|.+3282557.+|.+3294386.+|.+3294894.+|.+3303337.+|.+3320610.+|.+3335540.+|.+3344952.+|.+3345118.+|.+3355977.+|.+3401871.+|.+3443779.+|.+3534441.+|.+3591297.+|.+3595370.+|.+3601923.+|.+3624792.+|.+3689917.+|.+3689917.+|.+3695328.+|.+3696503.+|.+3762340.+|.+3775163.+|.+3870048.+|.+3924314.+|.+3935847.+|.+3940550.+',X1)) %>%
  write_delim(file = str_glue("{some_path}some_file_v1.vcf"),
                                 k                        col_names = F, quote = 'none', escape = 'none')

#zip the file
call = str_glue("bcftools view {some_path}some_file_v1.vcf -o {some_path}some_file_v1.vcf.gz")
system(call)
#index the new file
call = str_glue("bcftools index -f {f}.gz")
system(call)

Q1: How can one fix the large amount of REF prefix differences in the chromosome files?
Q2: If one uses the fixref plugin of bcftools which is the correct reference file to use after lifting was performed?
EDIT: JRodrigoF recommended to use the bcftools norm --check-ref s --fasta-ref GRCh38_full_analysis_set_plus_decoy_hla.fa  to fix the issue. It seems this produces the following error on my end:
 [E::faidx_adjust_position] The sequence "1" was not found faidx_fetch_seq failed at 1:69869
To make sure this is not a problem created during liftover I will post my liftover procedure here:
library(tidyverse)

#set working directory to data directory
trio_wd <- str_glue(here::here(),'/trio/SomeName/')

#create file list for raw data
file_list <- str_c(trio_wd,dir(trio_wd)) %>% str_extract('.+txt') %>%
  {.[!is.na(.)]}

    #generate map files from 23andme text file
    for(f in file_list){
      call = str_glue('plink --23file {f} --recode --out {str_remove(f,".txt")}
                      ')
      system(call)
    }

    #create file list for map data
    file_list <- str_c(trio_wd,dir(trio_wd)) %>% str_extract('.+map') %>%
      {.[!is.na(.)]}

    for (f in file_list){
    call <- str_glue('python3 reference/liftOverPlink.py --map {f} --out {str_replace(f,".map","_lifted")} --chain reference/hg19ToHg38.over.chain.gz  ')
    system(call)

    #create valid mapping file for hg19
    call <- str_glue('python3 reference/rmBadLifts.py -m {str_replace(f,".map","_lifted")}.map --out {str_replace(f,".map","_good_lifted")}.map -l {str_replace(f,".map","_bad_lifted.dat")}')
    system(call)

    call <- str_glue('cut -f 2 {str_replace(f,".map","_bad_lifted.dat")} > {str_replace(f,".map","_exclude_lifted.dat")}')
    system(call)
    call <- str_glue('cut -f 4 {str_replace(f,".map","_lifted")}.bed.unlifted | sed "/^#/d" >> {str_replace(f,".map","_exclude_lifted.dat")} ')
    system(call)
    # Note: this will clobber the lifted MAP file generated by `liftOverPlink`:
    call <- str_glue('plink --file {str_remove(f,".map")} --recode --out {str_replace(f,".map","_lifted")} --exclude {str_replace(f,".map","_exclude_lifted.dat")} ')
    system(call)
    call <- str_glue('plink --ped {str_replace(f,".map","_lifted.ped")} --map {str_replace(f,".map","_good_lifted")}.map --recode --out {str_replace(f,".map","_final")}')
    system(call)
    }

EDIT II:
Sample of the file including header
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##fileDate=20221108
##source=PLINKv1.90
##contig=<ID=1,length=248928329>
##contig=<ID=2,length=242157857>
##contig=<ID=3,length=198157342>
##contig=<ID=4,length=190106194>
##contig=<ID=5,length=181288140>
##contig=<ID=6,length=170595721>
##contig=<ID=7,length=159331484>
##contig=<ID=8,length=145075470>
##contig=<ID=9,length=138211490>
##contig=<ID=10,length=133659511>
##contig=<ID=11,length=135075227>
##contig=<ID=12,length=133261768>
##contig=<ID=13,length=114341522>
##contig=<ID=14,length=106873816>
##contig=<ID=15,length=101922277>
##contig=<ID=16,length=90172890>
##contig=<ID=17,length=83203771>
##contig=<ID=18,length=80252738>
##contig=<ID=19,length=58586567>
##contig=<ID=20,length=64323573>
##contig=<ID=21,length=46679699>
##contig=<ID=22,length=50776369>
##INFO=<ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on real reference genome">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##bcftools_viewVersion=1.13+htslib-1.13+ds
##bcftools_viewCommand=view -r 1 -o /home/oem/Documents/GitHub/Genomics_prac_guide/trio/some_name/some_name//genome_some_name_v5_Full_20191220113716_final_1.vcf.gz /home/oem/Documents/GitHub/Genomics_prac_guide/trio/some_name/some_name_v5_Full_20191220113716_final.vcf.gz; Date=Tue Nov  8 13:02:01 2022
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  FAM001_ID001
1   69869   rs548049170 T   .   .   .   PR  GT  0/0
1   74792   rs13328684  N   .   .   .   PR  GT  ./.
1   792461  rs116587930 G   .   .   .   PR  GT  0/0
1   817341  rs3131972   G   .   .   .   PR  GT  0/0
1   818725  rs12184325  C   .   .   .   PR  GT  0/0
1   820888  rs12567639  A   .   .   .   PR  GT  0/0
1   823656  rs114525117 G   .   .   .   PR  GT  0/0
1   858952  rs12127425  G   .   .   .   PR  GT  0/0
1   866156  rs79373928  T   .   .   .   PR  GT  0/0
1   880041  rs72888853  N   .   .   .   PR  GT  ./.
1   889018  rs7538305   A   .   .   .   PR  GT  0/0
1   894801  rs28444699  A   .   .   .   PR  GT  0/0
1   895351  i713449 N   .   .   .   PR  GT  ./.
1   899450  rs116452738 G   .   .   .   PR  GT  0/0
1   899712  rs72631887  T   .   .   .   PR  GT  0/0
1   903285  rs28678693  T   C   .   .   PR  GT  0/1
1   905373  rs4970382   T   C   .   .   PR  GT  0/1
1   911428  rs4475691   C   .   .   .   PR  GT  0/0
1   916010  rs72631889  G   .   .   .   PR  GT  0/0

EDIT III:
Output of grep '^>' reference/GRCh38_full_analysis_set_plus_decoy_hla.fa
>chr1  AC:CM000663.2  gi:568336023  LN:248956422  rl:Chromosome  M5:6aef897c3d6ff0c78aff06ac189178dd  AS:GRCh38
>chr2  AC:CM000664.2  gi:568336022  LN:242193529  rl:Chromosome  M5:f98db672eb0993dcfdabafe2a882905c  AS:GRCh38
>chr3  AC:CM000665.2  gi:568336021  LN:198295559  rl:Chromosome  M5:76635a41ea913a405ded820447d067b0  AS:GRCh38
>chr4  AC:CM000666.2  gi:568336020  LN:190214555  rl:Chromosome  M5:3210fecf1eb92d5489da4346b3fddc6e  AS:GRCh38
>chr5  AC:CM000667.2  gi:568336019  LN:181538259  rl:Chromosome  M5:a811b3dc9fe66af729dc0dddf7fa4f13  AS:GRCh38  hm:47309185-49591369

EDIT IV: Edited the fasta-file and the corresponding index file as mentioned in the comments. Still, the error from EDIT I sustains.
First five chromosomes of the modified FASTA-file:
grep '^>' reference/GRCh38_full_analysis_set_plus_decoy_hla_modified.fa
>1  AC:CM000663.2  gi:568336023  LN:248956422  rl:Chromosome  M5:6aef897c3d6ff0c78aff06ac189178dd  AS:GRCh38
>2  AC:CM000664.2  gi:568336022  LN:242193529  rl:Chromosome  M5:f98db672eb0993dcfdabafe2a882905c  AS:GRCh38
>3  AC:CM000665.2  gi:568336021  LN:198295559  rl:Chromosome  M5:76635a41ea913a405ded820447d067b0  AS:GRCh38
>4  AC:CM000666.2  gi:568336020  LN:190214555  rl:Chromosome  M5:3210fecf1eb92d5489da4346b3fddc6e  AS:GRCh38
>5  AC:CM000667.2  gi:568336019  LN:181538259  rl:Chromosome  M5:a811b3dc9fe66af729dc0dddf7fa4f13  AS:GRCh38  hm:47309185-49591369

Addjusted the bcftools-call of EDIT I afterwards so the modified FASTA-file will be used
call = str_glue('bcftools  norm --check-ref s {f} --fasta-ref reference/GRCh38_full_analysis_set_plus_decoy_hla_modified.fa -o {out} ')
[E::faidx_adjust_position] The sequence "1" was not found
faidx_fetch_seq failed at 1:69869

View region in question
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##fileDate=20221108
##source=PLINKv1.90
##contig=<ID=1,length=248928329>
##contig=<ID=2,length=242157857>
##contig=<ID=3,length=198157342>
##contig=<ID=4,length=190106194>
##contig=<ID=5,length=181288140>
##contig=<ID=6,length=170595721>
##contig=<ID=7,length=159331484>
##contig=<ID=8,length=145075470>
##contig=<ID=9,length=138211490>
##contig=<ID=10,length=133659511>
##contig=<ID=11,length=135075227>
##contig=<ID=12,length=133261768>
##contig=<ID=13,length=114341522>
##contig=<ID=14,length=106873816>
##contig=<ID=15,length=101922277>
##contig=<ID=16,length=90172890>
##contig=<ID=17,length=83203771>
##contig=<ID=18,length=80252738>
##contig=<ID=19,length=58586567>
##contig=<ID=20,length=64323573>
##contig=<ID=21,length=46679699>
##contig=<ID=22,length=50776369>
##INFO=<ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on real reference genome">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##bcftools_viewVersion=1.13+htslib-1.13+ds
##bcftools_viewCommand=view -r 1 -o /home/oem/Documents/GitHub/Genomics_prac_guide/trio/Some_name/Some_name//genome_Some_name_v5_Full_20191220113716_final_1.vcf.gz /home/oem/Documents/GitHub/Genomics_prac_guide/trio/Some_name/genome_Some_name_v5_Full_20191220113716_final.vcf.gz; Date=Tue Nov  8 13:02:01 2022
##bcftools_viewCommand=view -o /home/oem/Documents/GitHub/Genomics_prac_guide/trio/Some_name/Some_name/genome_Some_name_v5_Full_20191220113716_final_1_1.vcf.gz /home/oem/Documents/GitHub/Genomics_prac_guide/trio/Some_name/Some_name/genome_Some_name_v5_Full_20191220113716_final_1_1.vcf; Date=Tue Nov  8 18:16:55 2022
##bcftools_viewCommand=view -r 1:69869 /home/oem/Documents/GitHub/Genomics_prac_guide/trio/Some_name/Some_name/genome_Some_name_v5_Full_20191220113716_final_1_1.vcf.gz; Date=Thu Nov 10 10:56:54 2022
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  FAM001_ID001
1   69869   rs548049170 T   .   .   .   PR  GT  0/0

What about using bctools norm instead to fix the REF allele whenever needed? bcftools norm --check-ref s --fasta-ref GRCh38_full_analysis_set_plus_decoy_hla.fa
Or any other flavour of GRCh38 you are using, since the issue to solve is after the lifting over, so v38 has to be used in this case

If I do this I will recieve the following error msg: [E::faidx_adjust_position] The sequence "1" was not found faidx_fetch_seq failed at 1:69869 Maybe something went wrong during the lifting? I will add my lifting procedure to the OP.

looks to me that the header of your lifted VCF now in hg38 is not compatible with the sequence names contained in your reference fasta file; for example if in the hg38 version you have names like "chr1", ... and in the new VCF you have "1", an error like this would come up,

Added a sample of the file including the header. Does not seem that this is the case for me.

That error just means you have 1 in your VCF and probably chr1 in your fasta sequence.

How can I solve this error?

can you write here or include in the post the names of the chromosomes in your fasta file (reference). You can fetch the first 5 with grep '^>' <my_reference_file>

Done. So I guess I have to edit the fasta file and remove the "chr" string infront of each position?

Yes. You can do things like this with sed command. sed 's/>chr/>/' reference/GRCh38_full_analysis_set_plus_decoy_hla.fa > reference/GRCh38_full_analysis_set_plus_decoy_hla_modified.fa
Also I would check if the plink lifter does not have any relevant option regarding this issue. My expectation would be that the lifter produces a new VCF that contains the correct contig names as taken from the new reference, ... but here is not being the case. I would usually go for the lifter tool from GATK for example,

Done. Also adjusted the index file of the fast file in the same way. Still, I receive the same error as before. I added the the output of the region in question to the section EDIT IV of the question.

Good. I'll wrap it up as an answer too and maybe you could mark it as such so that I earn a few points,

It is important that after lifting over a VCF/BCF file from one coordinate system to another, the header of the resulting file contains the chromosome/contig names as present in the coordinate system of the new reference genome. it works as a double check that the assumptions of the coordinate system (particular reference) are consistent with the actual data (VCF lines). PLINK files for example do not have a header section and this doesn't not apply.
Looks like BCFTools (merge) complains whenever if finds positions that have one base ('C') as reference in one of the input files, and another base ('N') in another file. When working with VCF/BCF files, bcftools norm can be used to fix 'N' reference position into their actual base (if that is the case) and/or fix bases listed as reference (e.g. coming from a PLINK major-allele based BED file), but which actually correspond to an alternate allele based on either GRCh37/GRCh38, ... .
bcftools norm --check-ref s --fasta-ref <reference_file> <my_VCF_file>

Thanks for @JRodrigoF for pointing towards bcftools norm --check-ref s --fasta-ref. This will indeed fix the initial REF prefix diff issue of my question.
If anyone else stumbles over this question while facing a similar error to :[E::faidx_adjust_position] The sequence "1" was not found faidx_fetch_seq failed at 1:69869
Please try to re-index your FASTA-file using samtools
samtools faidx reference_file.fa Most of the times this seems to solve this kind of error!


text:
eQTL data extraction from database JSON file
I have a long list of JSON data, with repeats of contents similar to followings.
Due to the original JSON file is too long, I will just shared the hyperlinks here. This is a result generated from a database called RegulomeDB.
Direct link to the JSON file
I would like to extract specific data (eQTLs) from "method": "eQTLs" and "value": "xxxx", and put them into 2 columns (tab delimited) exactly like below.
Note: "value":"xxxx" is extracted right after "method": "eQTLs"is detected.
eQTLs   firstResult, secondResult, thirdResult, ...

In this example, the desired output is:
eQTLs   EIF3S8, EIF3CL

I've tried using a python script but was unsuccessful.
import json
with open('file.json') as f:
    f_json = json.load(f)
    print 'f_json[0]['"method": "eQTLs"'] + "\t" + f_json[0]["value"]

Thank you for your kind help.

Could you clarify what output format you expect please? You say "tab delimited" but show comma delimited data.

HI @terdon, edited my post. Thank you.

Here a jq solution if you can use it:
jq -r --arg f "eQTLs" '$f + "\t" + ([."@graph"[] | select(.method == $f) | .value] | join(", "))' file.json

Explanations:

-r: Raw string result (so result string will be without between ")
--arg name value: Your filter on eQTLs method value
Search in @graph array elements with []
Select only array element with method attribute equals to eQTLs (in $f argument)
The result, 2 in this case, are put in a array [...] and join with string ,
$f string, a tabulation "\t" and the join result are concatenated to one string with + operators

Result:
eQTLs   EIF3S8, EIF3CL

It's a little hard to know what your goal is here based on the sample code, but I do find that your JSON snippet does not match any spec. I was able to edit it to add brackets to make the entries elements of a list:
[
  {
  "biosample_ontology": {
    "@id": "https://www.encodeproject.org/biosample-types/primary_cell_CL_0000576/",
    "@type": [
      "BiosampleType",
      "Item"
    ],
    "aliases": [],
    "cell_slims": [
      "hematopoietic cell",
      "monocyte",
      "myeloid cell",
      "leukocyte"
    ],
    "classification": "primary cell",
    "dbxrefs": [],
    "developmental_slims": [
      "mesoderm"
    ],
    "name": "primary_cell_CL_0000576",
    "organ_slims": [
      "blood",
      "bodily fluid"
    ],
    "references": [],
    "schema_version": "1",
    "status": "released",
    "synonyms": [],
    "system_slims": [
      "immune system"
    ],
    "term_id": "CL:0000576",
    "term_name": "monocyte",
    "uuid": "20685264-1187-4c5c-b4ce-dc85ea53fc11"
  },
  "chrom": "chr11",
  "dataset": "https://www.encodeproject.org/annotations/ENCSR464XOM/",
  "dataset_rel": "/annotations/ENCSR464XOM/",
  "documents": [],
  "end": 59936979,
  "file": "ENCFF814NRI",
  "method": "eQTLs",
  "start": 59936978,
  "strand": null,
  "targets": [],
  "value": "MS4A4A"
},

    {
  "biosample_ontology": {
    "@id": "https://www.encodeproject.org/biosample-types/primary_cell_CL_0000576/",
    "@type": [
      "BiosampleType",
      "Item"
    ],
    "aliases": [],
    "cell_slims": [
      "hematopoietic cell",
      "monocyte",
      "myeloid cell",
      "leukocyte"
    ],
    "classification": "primary cell",
    "dbxrefs": [],
    "developmental_slims": [
      "mesoderm"
    ],
    "name": "primary_cell_CL_0000576",
    "organ_slims": [
      "blood",
      "bodily fluid"
    ],
    "references": [],
    "schema_version": "1",
    "status": "released",
    "synonyms": [],
    "system_slims": [
      "immune system"
    ],
    "term_id": "CL:0000576",
    "term_name": "monocyte",
    "uuid": "20685264-1187-4c5c-b4ce-dc85ea53fc11"
  },
  "chrom": "chr11",
  "dataset": "https://www.encodeproject.org/annotations/ENCSR464XOM/",
  "dataset_rel": "/annotations/ENCSR464XOM/",
  "documents": [],
  "end": 59936979,
  "file": "ENCFF814NRI",
  "method": "eQTLs",
  "start": 59936979,
  "strand": null,
  "targets": [],
  "value": "MS4A4B"
}
]

When I did this, I was able to slightly modify your code to print out each element's method key, value pair:
import json
with open('/Users/max.press/Downloads/file.json') as f:
    f_json = json.load(f)
    for i in range(len(f_json)):
        print(f_json[i]["method"] + "\t" + f_json[i]["value"])

% python ~/Downloads/readfile.py
# prints this:
eQTLs   MS4A4A
eQTLs   MS4A4B

Hi Maximilian Press, I am sorry the code doesn't work, with error code as followings.
print(f_json[i]["method"] + "\t" + f_json[i]["value"])
KeyError: 0

I am sorry about the JSON snippet format. Let me reupload a complete one to the original question. Could you please try again? Thank you.

@austin7923 it looks like there was additional formatting in your json outside of the snippet that you first posted, which makes the outer layer into a dict within which the list is nested. Thus the keyerror. regardless it looks like the jq solution works for you.

A (very) simple Python solution:
The script file readfile.py:
#! /usr/bin/python3

import sys
import json

filename=sys.argv[1]
filter=sys.argv[2]

results=[]
with open(filename) as json_file:
    json_content = json.load(json_file)
    for json_graph_element in json_content["@graph"]:
        if json_graph_element["method"] == filter:
            results.append(json_graph_element["value"])
print("%s\t%s" % (filter, ", ".join(results)))

Make the script executable:
chmod +x readfile.py

Use the script:
./readfile.py /path/to/the/file.json  eQTLs

Result:
eQTLs   EIF3S8, EIF3CL


text:
Multivariate analysis of cell culture bioprocess
I am seeking statistical support for the following bacterial growth time series.
In the lab, we have cell cultures growing in bioreactors. Every day we measure the concentration in mmol-1 of a set of metabolites and amino acids in the media with a cell culture analyser.
Some metabolites are consumed by the cells and some are produced during the time-course of the culture. Along with this, we also count the number of cells every day to calculate the cell growth dynamics.
We wish to assess the relationship between the consumption or production of any of these compounds and cell growth. For example, if the consumption of a specific metabolite from the media has a significant effect on increasing the cell count and would assist media formulation.
One measurement per day is taken over four consecutive days from multiple bioreactors.
What kind of analysis is would be appropriate?

This is a machine learning application. Its clear cut. Whether you want to perform unsupervised learning is your call.
Input data is for the cell growth, metabolite concentration time series would be:

T1 ... T4 (cell growth) + T1 ... T4 (metabolite 1) etc ...

The training target is:

media type

I would start with a simple naive Bayes model and see how it goes. If you obtain a high accuracy score you're in business, its worked. You would simply plot the result. The classification may not be linear in which case an alternative model is used.
You can manually test it, by inputing the data of cell growth etc ... and seeing if it correctly predicts the media type.
You then ask the ML algorithm which the metabolite/cell growth the 'media type' is influencing and it will tell you. Setting it up is complex I agree, but thats the gist and was the question.
Just to finish, if you want to limit what a given media is in particular influencing, i.e. not loads of measurements but just a few inputs/variables, e.g. a metabolite at a particular time point and/or cell growth at a given time point, then ML delivers that information. Thats a key part of its strength. Thus you can then trial loads of different media - even stuff not trained - and just measure a small number of select inputs/variables rather than the entire kitchen sink - which is what the OP is doing right now.

The overall input of ML into biological problems is a bit worrying in truth and tends to be over-represented in a select number of fields, protein modelling, and zero in everything else.


text:
Compare and Reorganize Fasta Headers Python
I want to compare the headers from the fasta file to file1, and if there's a match, reorganize the header and put the match first. If there's no match between fasta file and file1,  look at file2 and if there's a match, put that name first on the header.
The goal is to organize the fasta file header with the similar match listed as the first entry on the header. How can I do this with this code or is there another way to write it in python or biopython?
BioPython solutions welcome.
Current files look like this
seqs.fasta:
>ant| bee| fly
AGCT...
>bird| eagle| hawk| vulture
GATCG...

file1.txt:
mouse
hawk
cow

file2.txt:
crane
fly

The output file should look like this, so that the order of the header starts with the match from file 1 or file 2:
> fly| ant| bee
AGCT...
> hawk| bird| eagle| vulture
AGCT...

My code is:
fasta= open('seq.fasta', 'r')
one= open('file1.txt' , 'r')
two= open('file2.txt' , 'r')
output= open('newseq.fasta', 'w')

for line in fasta:
    if line.startswith('>'): # getting just the header
        header= line.strip().split("|") # there are several entries on the header separated by the |
        for name in one:
           if name in header: # looking to see if file1.txt has match in fasta file
              header.insert(0,line.pop(header.index(name))) # insert the name from file1 as the first header entry
        else: # if no match in file1 look at file2
            for name in two:
               if name in header:
                    header.insert(0,line.pop(header.index(name))) # insert the name from file2 as first header entry
    else:
output.write(line)

There is nothing wrong with the concept of the code, that is exactly how I would do it. The problem is here:
    else:
output.write(line)

Corrections
import os
output = 'newseq.fasta'
if os.path.exists(output):
    os.unlink(output)

for line in fasta:
    if line.startswith('>'): # getting just the header
        header= line.strip().split("|") # there are several entries on the header separated by the |
        for name in one:
           if name in header:
              line = header.insert(0,line.pop(header.index(name))) # insert the name from file1 as the first header entry
        else: # if no match in file1 look at file2
            for name in two:
               if name in header:
                    line = header.insert(0,line.pop(header.index(name)))

# The indentation must be absolutely correct
# otherwise the code falls over. Must be one indent from
# the first for statement
    with open(output, 'a') as fout:
        fout.write(line)

I assume that fasta, one and two are loaded successfully.

Thanks :-) @NoraJob upvoted previously.


text:
Conversion FASTA to BED fails
I want to use conform-gt to get rid of some REF prefix differences see here. Conform-gt needs the reference FASTA file split by chromosome and in .vcf.gz format. How can I do this?
I already tried faidx with
faidx --transform bed reference/GRCh38_full_analysis_set_plus_decoy_hla.fa > reference/GRCh38_full_analysis_set_plus_decoy_hla.bed

but this produced the following error:
Traceback (most recent call last):
  File "/home/oem/.local/bin/faidx", line 8, in <module>
    sys.exit(main())
  File "/home/oem/.local/lib/python3.9/site-packages/pyfaidx/cli.py", line 202, in main
    write_sequence(args)
  File "/home/oem/.local/lib/python3.9/site-packages/pyfaidx/cli.py", line 51, in write_sequence
    outfile.write(transform_sequence(args, fasta, name, start, end))
  File "/home/oem/.local/lib/python3.9/site-packages/pyfaidx/cli.py", line 121, in transform_sequence
    line_len = fasta.faidx.index[name].lenc
KeyError: 'HLA-A*01'

Can you please [edit] your question and explain what you mean a bit more? Your title talks about converting fasta (sequences) to bed (regions) which is odd, but could be done if you mean "read the fasta and tell me what regions of a reference genome they map to". But then your question talks about converting fasta to "vcf.gz" format and this doesn't make sense unless you want some sort of gVCF file with an entry for every position of a genome.

I reconsidered my question. I don't think that I really fully understood what I want to achieve with this format-conversion nor if I really understood what conform-gt really does. Back to the drawing board I guess. I will close this question and hope that the original question in my OP can solve my problem.


text:
finding matching motifs in a sequence and their position
I am trying to find some matching motifs on a sequence, as well as the position that the motif is located in and then output that into a fasta file. The code below shows that the motif [L*L*L] is present in the sequence, when I run it returns as "YES" but I do not know where it is positioned.
peptide1="MKFSNEVVHKSMNITEDCSALTGALLKYSTDKSNMNFETLYRDAAVESPQHEVSNESGSTLKEHDYFGLSEVSSSNSSSGKQPEKCCREELNLNESATTLQLGPPAAVKPSGHADGADAHDEGAGPENPAKRPAHHMQQESLADGRKAAAEMGSFKIQRKNILEEFRAMKAQAHMTKSPKPVHTMQHNMHASFSGAQMAFGGAKNNGVKRVFSEAVGGNHIAASGVGVGVREGNDDVSRCEEMNGTEQLDLKVHLPKGMGMARMAPVSGGQNGSAWRNLSFDNMQGPLNPFFRKSLVSKMPVPDGGDSSANASNDCANRKGMVASPSVQPPPAQNQTVGWPPVKNFNKMNTPAPPASTPARACPSVQRKGASTSSSGNLVKIYMDGVPFGRKVDLKTNDSYDKLYSMLEDMFQQYISGQYCGGRSSSSGESHWVASSRKLNFLEGSEYVLIYEDHEGDSMLVGDVPWELFVNAVKRLRIMKGSEQVNLAPKNADPTKVQVAVG"

if re.search(r"L*L*L", peptide1):
    print("YES")
else:
    print("NO")

The code that I wrote to find the position is below, but when I run it says invalid syntax. Could you please assist as I am still new in Bioinformatics, and I have no clue whether in the right track or not.
for position in range(len(peptide1)):
    if peptide1[position:].startswith(re.search(r"L*L*L")):
        print(position+1)

Lastly once the position of these motifs has been identified, how do I highlight them such that they are visible in the output fasta file.
I am searching for LxLxL where x can be any amino acid. I used the * to indicate to python that it should use any amino acid, but I think I am incorrect. Could you please assist on how I should actually be writing that motif?

What is the actual pattern you are looking for? L*L*L* will match anything at all since L* means "0 or more L", so it will match L but it will also match X or anything else.

@terdon, I am searching for LxLxL where x can be any amino acid. I used the * to indicate to python that it should use any amino acid , but I think I am incorrect. Could you please assist on how I should actually be writing that motif.

@terdon If for example I have a motif LLxxxLxPP, could you please advice on how I am supposed to write it inside re.search or re.finditer? this is where I get confused

re.search method do not return a simply boolean but a useful object contains many information:
#! /usr/bin/env python3

import re

peptide1="MKFSNEVVHKSMNITEDCSALTGALLKYSTDKSNMNFETLYRDAAVESPQHEVSNESGSTLKEHDYFGLSEVSSSNSSSGKQPEKCCREELNLNESATTLQLGPPAAVKPSGHADGADAHDEGAGPENPAKRPAHHMQQESLADGRKAAAEMGSFKIQRKNILEEFRAMKAQAHMTKSPKPVHTMQHNMHASFSGAQMAFGGAKNNGVKRVFSEAVGGNHIAASGVGVGVREGNDDVSRCEEMNGTEQLDLKVHLPKGMGMARMAPVSGGQNGSAWRNLSFDNMQGPLNPFFRKSLVSKMPVPDGGDSSANASNDCANRKGMVASPSVQPPPAQNQTVGWPPVKNFNKMNTPAPPASTPARACPSVQRKGASTSSSGNLVKIYMDGVPFGRKVDLKTNDSYDKLYSMLEDMFQQYISGQYCGGRSSSSGESHWVASSRKLNFLEGSEYVLIYEDHEGDSMLVGDVPWELFVNAVKRLRIMKGSEQVNLAPKNADPTKVQVAVG"

re_search_result = re.search(r'L*L*L', peptide1)
if re_search_result:
    print("YES")
    motif_start_position=re_search_result.start()
    motif_end_position=re_search_result.end()
    motif_length = re_search_result.end() - re_search_result.start()
    print("Motif start at position: %d" % motif_start_position)
    print("Motif end at position  : %d" % motif_end_position)
    print("SO length is ......... : %d" % motif_length)
    print("The motif is : <%s>" % peptide1[motif_start_position:motif_end_position])
else:
    print("NO")

Result is:
YES
Motif start at position: 20
Motif end at position  : 488
SO length is ......... : 468
The motif is : <LTGALLKYSTDKSNMNFETLYRDAAVESPQHEVSNESGSTLKEHDYFGLSEVSSSNSSSGKQPEKCCREELNLNESATTLQLGPPAAVKPSGHADGADAHDEGAGPENPAKRPAHHMQQESLADGRKAAAEMGSFKIQRKNILEEFRAMKAQAHMTKSPKPVHTMQHNMHASFSGAQMAFGGAKNNGVKRVFSEAVGGNHIAASGVGVGVREGNDDVSRCEEMNGTEQLDLKVHLPKGMGMARMAPVSGGQNGSAWRNLSFDNMQGPLNPFFRKSLVSKMPVPDGGDSSANASNDCANRKGMVASPSVQPPPAQNQTVGWPPVKNFNKMNTPAPPASTPARACPSVQRKGASTSSSGNLVKIYMDGVPFGRKVDLKTNDSYDKLYSMLEDMFQQYISGQYCGGRSSSSGESHWVASSRKLNFLEGSEYVLIYEDHEGDSMLVGDVPWELFVNAVKRLRIMKGSEQVNL>

The motif L*L*L is a quite strange: L zero or many times followed by L zero or many times followed by followed by a L. It's an equivalent of L*L, equivalent of L+.
Do you search for L.*L.*L, L[^L]*L[^L]*L or simply LLL ?
Notes:

When you post code, use block delimited by two ``` (no by two simple `)
In Python, code indentation is very important, always use spaces (4)

If you want to loop over all motifs (like strings), use findall method:
#! /usr/bin/env python3

import re

peptide1="MKFSNEVVHKSMNITEDCSALTGALLKYSTDKSNMNFETLYRDAAVESPQHEVSNESGSTLKEHDYFGLSEVSSSNSSSGKQPEKCCREELNLNESATTLQLGPPAAVKPSGHADGADAHDEGAGPENPAKRPAHHMQQESLADGRKAAAEMGSFKIQRKNILEEFRAMKAQAHMTKSPKPVHTMQHNMHASFSGAQMAFGGAKNNGVKRVFSEAVGGNHIAASGVGVGVREGNDDVSRCEEMNGTEQLDLKVHLPKGMGMARMAPVSGGQNGSAWRNLSFDNMQGPLNPFFRKSLVSKMPVPDGGDSSANASNDCANRKGMVASPSVQPPPAQNQTVGWPPVKNFNKMNTPAPPASTPARACPSVQRKGASTSSSGNLVKIYMDGVPFGRKVDLKTNDSYDKLYSMLEDMFQQYISGQYCGGRSSSSGESHWVASSRKLNFLEGSEYVLIYEDHEGDSMLVGDVPWELFVNAVKRLRIMKGSEQVNLAPKNADPTKVQVAVG"

re_findall_result = re.findall(r'L[^L]*L[^L]*L', peptide1)
for element in re_findall_result:
    print("Element: <%s>" % element)

Result will be:
Element: <LTGALL>
Element: <LYRDAAVESPQHEVSNESGSTLKEHDYFGL>
Element: <LNLNESATTL>
Element: <LGPPAAVKPSGHADGADAHDEGAGPENPAKRPAHHMQQESLADGRKAAAEMGSFKIQRKNIL>
Element: <LDLKVHL>
Element: <LSFDNMQGPLNPFFRKSL>
Element: <LVKIYMDGVPFGRKVDLKTNDSYDKL>
Element: <LEDMFQQYISGQYCGGRSSSSGESHWVASSRKLNFL>
Element: <LIYEDHEGDSMLVGDVPWEL>

If you want all find objects like with search method, use finditer method:
#! /usr/bin/env python3

import re

peptide1="MKFSNEVVHKSMNITEDCSALTGALLKYSTDKSNMNFETLYRDAAVESPQHEVSNESGSTLKEHDYFGLSEVSSSNSSSGKQPEKCCREELNLNESATTLQLGPPAAVKPSGHADGADAHDEGAGPENPAKRPAHHMQQESLADGRKAAAEMGSFKIQRKNILEEFRAMKAQAHMTKSPKPVHTMQHNMHASFSGAQMAFGGAKNNGVKRVFSEAVGGNHIAASGVGVGVREGNDDVSRCEEMNGTEQLDLKVHLPKGMGMARMAPVSGGQNGSAWRNLSFDNMQGPLNPFFRKSLVSKMPVPDGGDSSANASNDCANRKGMVASPSVQPPPAQNQTVGWPPVKNFNKMNTPAPPASTPARACPSVQRKGASTSSSGNLVKIYMDGVPFGRKVDLKTNDSYDKLYSMLEDMFQQYISGQYCGGRSSSSGESHWVASSRKLNFLEGSEYVLIYEDHEGDSMLVGDVPWELFVNAVKRLRIMKGSEQVNLAPKNADPTKVQVAVG"

re_findall_result = re.finditer(r'L[^L]*L[^L]*L', peptide1)
for element in re_findall_result:
    print("peptide1[%s:%s]=<%s>" % (element.start(), element.end(), peptide1[element.start():element.end()]))

Result will be:
peptide1[20:26]=<LTGALL>
peptide1[39:69]=<LYRDAAVESPQHEVSNESGSTLKEHDYFGL>
peptide1[90:100]=<LNLNESATTL>
peptide1[101:163]=<LGPPAAVKPSGHADGADAHDEGAGPENPAKRPAHHMQQESLADGRKAAAEMGSFKIQRKNIL>
peptide1[248:255]=<LDLKVHL>
peptide1[278:296]=<LSFDNMQGPLNPFFRKSL>
peptide1[378:404]=<LVKIYMDGVPFGRKVDLKTNDSYDKL>
peptide1[407:443]=<LEDMFQQYISGQYCGGRSSSSGESHWVASSRKLNFL>
peptide1[449:469]=<LIYEDHEGDSMLVGDVPWEL>

thank you so much.  I have another question, if for example the motif is written as LLxxLxPP, where x can be any amino acid, how must I write it? I think that is one of the things that confuses me

So if I have LxxLxPP where x can be any amino acid, how must I write it inside the re.finditer ? This is where I get confused

@thole in regular expressions, * means "0 or more of the previous character". To match any character, you want .. Also, [ABC] means "find either A or B or C". So the actual pattern you want is r'L.L.L' but of course that will also match LLLLL.

@terdon thank you, that makes a lot of sense. So if I had a motif LxxLxP then it would be r'L..L.P'

Yes, exactly right @thole.

if you are interested by Regular Expression (RegEx) in Python, see https://www.w3schools.com/python/python_regex.asp


text:
Is there any book or reference that explains evolutionary dynamics and fixation probability?
I am a master's student in computational and mathematical biology, with a bachelor's degree in pure mathematics.
I am really struggling with my computational biology course, we are supposed to study problems related to evolutionary dynamics like:
1)the rock-paper-scissors game which explains the evolution of male mating strategies in side-blotched lizards
2)fixation probability of a new mutation.
Could someone link me to some books or references that might be helpful in solving such problems?

https://drive.google.com/file/d/1qdSrAeTNqJRSFYiU335EK4oj2mSv97Id/view?usp=sharing

Maynard Smith, J. (1989) Evolutionary Genetics. Oxford: Oxford University Press. ISBN 0-19-850231-1
Maynard Smith, J. and Harper, D. (2003) Animal Signals. Oxford University Press. ISBN 0-19-852685-7[16]

Then perform a PubMed search for Adam Eyre-Walker.

I'm fond of Joe Felsenstein's book on evolutionary genetics as a relatively recent source that presents a lot of these problems at a very simple level for the newer student, and then adds complexity gradually. Page 260 seems relevant.


text:
Is current AI able to generate RNA sequences of viruses?
After seeing a video of a combination of CRISPR and AI, and an article of someone who made two babies immune against HIV and still healthy, I wondered about something. If such complicated thing, to edit someone's genome, so that the person isn't likely to get other diseases or complications, while having a selected function, is (almost) possible, will it be possible currently to let AI generate RNA of virus particles like from a Generative Adversarial Network?

We show that these
tools capture important structures of the data and, when applied to designing probes
for protein binding microarrays (PBMs), allow us to generate new sequences whose
properties are estimated to be superior to those found in the training data.

https://arxiv.org/abs/1712.06148 In this research article, it is stated that the researchers were successful in training the model so that it would generate new sequences, while having a specific functioning protein. Note that this was from the time when ProGAN wasn't even invented yet, so they used a Wasserstein-GAN.

While we found these techniques to generally be well-behaved, we
did notice some undesirable cases during our explorations.

They concluded that even one mistake in nucleotides could mess the protein up. But now in 2022, we have StyleGAN-Ada (https://github.com/NVlabs/stylegan2-ada), which is far better in generating stuff, we've already (almost) overcome that stuff? So if you basically train such model on RNA of virus particles (most viruses have really short genome compared to an eukaryote cell, so a supercomputer is not required for training such thing, and there are databases online: Publicly available genome sequence database for viruses?), such as every mutation of SARS-CoV-2, HIV, ebola and so on, it would basically be possible to let it generate something that is as contageous as SARS-CoV-2, whereas as deadly as HIV? Keep in mind that one single generated sequence that functions, would already be devastating. It doesn't need to be 100% accurate, because after all we can just experiment with it and find out whether it is working or not as supposed to be.
It doesn't have to be a GAN:

We present three approaches: creating
synthetic DNA sequences using a generative adversarial network (GAN); a DNAbased variant of the activation maximization (“deep dream”) design method; and a
joint procedure which combines these two approaches together.

The RNA-sequence can be brought into a cell, and the cell will automatically create the components of the virus, which will become the virus ultimately.
Anyone can recall that this is possible? I have done my searches on Google, but I haven't found any article particularly describing this concept.
I hope this is on-topic, I don't know how the forums work. Thanks in advance.

I'm closing this question for multiple reasons, including: opinion based, speculation about dangerous experiments, doesn't restrict to a specific problem. The user name "Random user" doesn't give me much incentive to be lenient on this.

Reinforced learning is an interesting technology, but they are prone to extreme results - which are not necessarily at all biological feasible. Google can use it to outperform, or very possibly outcompete, the world's best chess player, because all of the rules are known. In biology we simply do not know all the rules.
To reiterate, the problem with reinforced learning in biology is not only do we not know all the rules, in fact in reality we know very few of the rules. HIV - why does it cause persistent infection? SARS-CoV-2 why doesn't the vaccine confer complete protection? Moreover, the GAN originator needs to know all of the biological known rules and apply them within the GAN and there are a lot of rules. Does the originator know the conformational change of env or spike and on which amino acid that are invariant? I doubt it. That information would be pivotal to generate AI sequence for a biologically feasible virus.
Even in the original Atari video game experiment * whilst GANS outcompeted humans by >1000%, equally they under-performed humans considerable for around >20% games tested. Thus even if all the rules are known GANs can mess-up a large proportion of the time.
These are simply hypotheses in an emerging technology within biology that require biological verifiable results and those proofs will be difficult, i.e. time consuming. This is not the case within video games, or board games.
In the first instance, GAN output can be investigated for biological validity against strongly established bioinformatics analysis - the stuff we use to assess whether there are errors in the data. However, there needs to be a lot incentive to perform this research at this point in the technological cycle. When the authors have something concrete, we might then reinvestigate the fundamental issues of their predictions.
*, the old Atari video games were subject to GANs and the computer then compared its performance in the video game against an average human. The results were highly variable, moreover humans are highly variable - which to my knowledge was ever investigated - so it was never at all clear the "the average human" represented in these tests.

A GAN isn't 'trained against [biological entity]' in the technical way supervised learning is trained. A GAN 'trains' against itself via a 'win/lose' quantification. Establishing a valid 'win/lose' is a pivotal part of the GAN. Thus in chess thats easy ... you win or lose the game, in biological that is not trivial.
Note
To address the point below, yes GANs devise their own "rules" to maximise win/loss WITHIN the rules already in place. However, they needs 'ground rules' otherwise how could it play a game of chess if it doesn't follow the rules of chess? How could it play an Atari video game if it does not know what a "joystick" is? The problem is without knowing the ground rules (what is biologically permitted), what's the validity of the subsequent rule set, i.e. predicted sequence in this case?
So a virus that can't attach to the host cellular receptor, isn't a virus because it will never be able to enter the cell.
Further note
The OP pointed out in the notes

"The idea of style transfer has largely only been explored in image-based tasks, which we attribute in part to the specific nature of loss functions used for style transfer."

The author is talking about a variant of transfer learning. What the author is reminding us is that overall transferring-learning is the technique which powered the success of computer vision. They are saying they have replicated that in their GANs.
Please note any further points should be separate questions

https://uzman-anwar.github.io/projects/IL_Atari/report.pdf Do you mean this article for atari?


text:
Turning GMT file into data frame so that third column is comma seperated
Similar to: How to convert data in gmt format to dataframe?
I have downloaded a GMT file from the GSEA database (https://www.gsea-msigdb.org/gsea/msigdb/download_file.jsp?filePath=/msigdb/release/2022.1.Hs/c2.all.v2022.1.Hs.symbols.gmt)
I have imported this into R using the GSA.read.gmt tool from GSA as I want to do pathway analysis using GAUSS (https://github.com/diptavo/GAUSS/blob/master/example_gmt.txt).
GAUSS needs the table to have the genes in a comma separated manner as well as there being headers in the table.
I have tried using the gsub function on the third column but it doesn't work, I feel this is because a GMT file is not the same as a dataframe.
Would anyone be able to help me turn this column into a comma separated one please?
Many thanks

I used the unite function from tidyr but first I read the gmt file in as a data table using the read.delim function (with header set to False). Then I used unite:
library(tidyr)

gmt_file <- read.delim("~/path_to_file/file.gmt", header = F)
gmt_file <- gmt_file %>% unite(genes,3:ncol(gmt_file), sep=",")
colnames(gmt_file) <- c("pathway", "source", "genes")
write.table....


text:
Converting Pairwise single- linkage clustering distance data to "newick" format
I have pairwise distance data for single linkage cluster and I would like to convert it to newick format.
I was unable to find a conversion algorithm.
I require newick format for visualization and annotation of a phylogeny.
Here is the output.
INFO  : ======================================
INFO  : Pairwise single linkage clustering
INFO  : ======================================
INFO  : Hierarchical Tree
INFO  : ======================================
INFO  : Node     Item 1   Item 2      Distance
INFO  :     0 :       53       52        0.001  5zxf  5z02
INFO  :    -1 :       64       60        0.002  7bj6  7bir
INFO  :    -2 :       62       -1        0.002  7biv
INFO  :    -3 :       34       26        0.002  4wt2  4ode
INFO  :    -4 :       65       -2        0.003  7bmg
INFO  :    -5 :       25       13        0.003  4occ  4erf
INFO  :    -6 :       29       -3        0.009  4ogt
INFO  :    -7 :       39       38        0.010  4zyi  4zyf
INFO  :    -8 :       50       48        0.011  5oc8  5ln2
INFO  :    -9 :       54       -8        0.012  6ggn
INFO  :   -10 :       42       41        0.015  5hmi  5hmh
INFO  :   -11 :       -4       61        0.020            7bit
INFO  :   -12 :       70       40        0.028  7qdq  5c5a
INFO  :   -13 :       47       68        0.032  5laz  7na3
INFO  :   -14 :       -9       45        0.035            5law
INFO  :   -15 :       49       59        0.038  5oai  6q9o
INFO  :   -16 :       32       -5        0.039  4qo4
INFO  :   -17 :      -15        8        0.044            3uu1
INFO  :   -18 :       23      -16        0.047  4oas
INFO  :   -19 :      -14       -7        0.053
INFO  :   -20 :       24      -18        0.057  4oba
INFO  :   -21 :      -19        6        0.057            3lbl
INFO  :   -22 :       43      -10        0.058  5hmk
INFO  :   -23 :       44       19        0.059  5lav  4jvr
INFO  :   -24 :      -23      -20        0.061
INFO  :   -25 :      -11       15        0.061            4hg7
INFO  :   -26 :       57      -24        0.062  6q9h
INFO  :   -27 :       27       -6        0.065  4odf
INFO  :   -28 :      -25      -17        0.067
INFO  :   -29 :       14      -28        0.067  4hbm
INFO  :   -30 :      -13      -21        0.067
INFO  :   -31 :      -27      -29        0.067
INFO  :   -32 :      -26      -30        0.068
INFO  :   -33 :       58       21        0.068  6q9l  4mdn
INFO  :   -34 :      -31      -32        0.070
INFO  :   -35 :       17        5        0.071  4jv9  3lbk
INFO  :   -36 :       30      -34        0.074  4ogv
INFO  :   -37 :      -33        2        0.077            1t4e
INFO  :   -38 :       35      -36        0.079  4zfi
INFO  :   -39 :       20      -38        0.086  4jwr
INFO  :   -40 :      -39      -37        0.094
INFO  :   -41 :        7      -40        0.106  3tj2
INFO  :   -42 :       51      -41        0.111  5trf
INFO  :   -43 :       69      -42        0.112  7na4
INFO  :   -44 :       56      -43        0.113  6q96
INFO  :   -45 :        4      -44        0.127  3jzk
INFO  :   -46 :       33      -45        0.129  4qoc
INFO  :   -47 :       22      -46        0.137  4mdq
INFO  :   -48 :      -12      -47        0.151
INFO  :   -49 :       37      -48        0.155  4zyc
INFO  :   -50 :       31      -49        0.156  4oq3
INFO  :   -51 :       12      -50        0.158  4ere
INFO  :   -52 :      -22       46        0.161            5lay
INFO  :   -53 :      -35      -51        0.182
INFO  :   -54 :       55      -53        0.186  6i29
INFO  :   -55 :       10      -54        0.224  3w69
INFO  :   -56 :       18      -55        0.238  4jve
INFO  :   -57 :        9      -56        0.265  3vzv
INFO  :   -58 :       28      -57        0.271  4ogn
INFO  :   -59 :      -58        1        1.000            1rv1
INFO  :   -60 :       66      -59        1.000  7na1
INFO  :   -61 :        0      -60        1.000
INFO  :   -62 :      -52      -61        1.000
INFO  :   -63 :       11      -62        1.000  4dij
INFO  :   -64 :       67      -63        1.000  7na2
INFO  :   -65 :       16      -64        1.000  4jv7
INFO  :   -66 :       36      -65        1.000  4zgk
INFO  :   -67 :        3      -66        1.000  2lzg

The problem with the tree format presented is that it does not conform to a known standard. Whilst its possible to write a parser I don't think its a good idea.
This output is hierarchical clustering, UPGMA, that sort of thing. Its not phylogenetics because this is a linkage tree.
Can you access the distance matrix prior to the clustering? If so the information below the line will help.
The output should provide the raw distance matrix prior to hierarchical clustering and will certainly have generated it in the calculation. BTW that particular clustering method isn't cool compared with nj (below) for complex reasons.
The calculation will be:

Perform pairwise linkage analysis
Perform cluster analysis on the resulting matrix

You need neighbor-joining to convert the distance matrix into a newick-format. There is no possibility of a direct conversion.
If you are using R its in the ape library via
treenk <- nj(M) with M is the matrix. The treenk will be a format that can be read straight into ggtree
Briefly, a distance matrix requires a clustering algorithm and neighbor-joining is by far the best clustering method developed in phylogenetics.
The following is a phylogenetics example ... this is different from the example you have 'cause thats about linkage.
data(mynucleotides)
tr <- nj(dist.dna(mynucleotides))
plot(tr)

I your case the final line will be something like ggtree(tr) ... with various options.

Thank for the idea. It helps so much.


text:
finding motifs in fasta file with multiple sequences
I am able to find matching motifs on a single sequence, as well as their position as shown below :
demo = "ABCLXLYLZLABC"
motif="L.L.L"
for element in re.finditer(motif,demo):
    print(element.start(),element.end,element.group())

This works and gives me the desired output which is:
 3 8 LXLYL
Now I want to know how can I then do this same process but on a multifasta file.
As done in this example script from https://www.biostars.org/p/14305/ .
The protein sequences on the file are :
> P1
MPPRRSIVEVKVLDVQKRRVPNKHYVYIIRVTWSSGATEAIYRRYSKFFDLQMQMLDKFP
MEGGQKDPKQRIIPFLPGKILFRRSHIRDVAVKRLIPIDEYCKALIQLPPYISQCDEVLQ
FFETRPEDLNPPKEEHIGKKKSGNDPTSVDPMVLEQYVVVADYQKQESSEISLSVGQVVD
IIEKNESGWWFVSTAEEQGWVPATCLEGQDGVQDEFSLQPEEEEKYTVIYPYTARDQDEM
NLERGAVVEVVQKNLEGWWKIRYQGKEGWAPASYLKKNSGEPLPPKLGPSSPAHSGALDL
DGVSRHQNAMGREKELLNNQRDGRFEGRLVPDGDVKQRSPKMRQRPPPRRDMTIPRGLNL
>P2
 MAEVRKFTKRLSKPGTAAELRQSVSEAVRGSVVLEKAKLVEPLDYENVITQRKTQIYSDP
 LRDLLMFPMEDISISVIGRQRRTVQSTVPEDAEKRAQSLFVKECIKTYSTDWHVVNYKYE
 DFSGDFRMLPCKSLRPEKIPNHVFEIDEDCEKDEDSSSLCSQKGGVIKQGWLHKANVNST
 ITVTMKVFKRRYFYLTQLPDGSYILNSYKDEKNSKESKGCIYLDACIDVVQCPKMRRHAF
 ELKMLDKYSHYLAAETEQEMEEWLIMLKKIIQINTDSLVQEKKDTVEAIQEEETSSQGKA
 ENIMASLERSMHPELMKYGRETEQLNKLSRGDGRQNLFSFDSEVQRLDFSGIEPDVKPFE
 EKCNKRFMVNCHDLTFNILGHIGDNAKGPPTNVEPFFINLALFDVKNNCKISADFHVDLN
 PPSVREMLWGTSTQLSNDGNAKGFSPESLIHGIAESQLCYIKQGIFSVTNPHPEIFLVVR
>P3
GDDSEWLKLPVDQKCEHKLWKARLSGYEEALKIFQKIKDEKSPEWSKYLGLIKKFVTDS
NAVVQLKGLEAALVYVENAHVAGKTTGEVVSGVVSKVFNQPKAKAKELGIEICLMYVEIE
KGESVQEELLKGLDNKNPKIIVACIETLRKALSEFGSKIISLKPIIKVLPKLFESRDKAV
RDEAKLFAIEIYRWNRDAVKHTLQNINSVQLKELEEEWVKLPTGAPKPSRFLRSQQELEA
KLEQQQSAGGDAEGGGDDGDEVPQVDAYELLDAVEILSKLPKDFYDKIEAKKWQERKEAL
EAVEVLVKNPKLEAGDYADLVKALKKVVGKDTNVMLVALAAKCLTGLAVGLRKKFGQYAG
HVVPTILEKFKEKKPQVVQALQEAIDAIFLTTTLQNISEDVLAVMDNKNPTIKQQTSLFI
ARSFRHCTSSTLPKSLLKPFCAALLKHINDSAPEVRDAAFEALGTALKVVGEKSVNPFLA

The script that was used to find the motif is below, the problem now is that it does not show the position where the motif is found. Could you assist in how one could do so.
import re
infile=open("seq.fasta",'r')
out=open("results.csv",'w')
pattern=re.compile(r"(P[A-Z]{2}P)")
for line in infile:
  line = line.strip("\n")
  if line.startswith('>'):
     name=line
else:
   s = re.findall(pattern,line)
   print '%s:%s' %(name,s)
   out.write('%s:\t%s\n' %(name,s))

Their script seems to have worked because the desired output was
>p1:
PGCP

>p1:
PHCP,
PKCP

The problem with their script is that, it does not provide the position where the sequence was found. Could you please assist on how that script can be modified? I want to try it but I want my output to have (1) Name Of Sequence, (2) The Motif Found , (3) The Position

This is easy just replace findall with finditer in the above code. The index will then appear where the match first appears and where it stops. However, you will to parse the finditer object which in the code is s
Then replace
print '%s:%s' %(name,s)

with
for match_obj in s:
    # print each re.Match object
    print(match_obj)

    # extract each matching number
    print(match_obj.group())

That should work. When it does you simply pass the match_obj and match_obj.group() to a file by replacing s with match_obj, match_obj.group()
Thus replace
  out.write('%s:\t%s\n' %(name,s))

with
  out.write('%s:\t%s\t%s\n' %(name,match_obj,match_obj.group()))

You might need to do an append BTW.
thus replace
out=open("results.csv",'w')

with
import os
if os.path.exists("results.csv"):
    os.unlink ("results.csv")
out=open("results.csv",'a')

Thank you so much for your suggestion, I am still new python user and will read up on biopython. With regards to the code will it now be: else:s=re.finditer (pattern,line)
print'%s:%' %(name,s)  ? Also by parse the finditer object what does that mean? I am struggling to link it to the original code. Could you perhaps try to edit your answer such that it looks like the original code, so that I can see where the changes will need to be made.


text:
Calling isoforms from long read data generated from partially degraded RNA
What will be the best tool to call isoforms from long read data generated from partially degraded RNA. By mistake we processed some samples with poor quality RNA to generate long read. Now we are looking for some tool which can merged the reads for the potential degraded transcripts. I think majority of the available tools like Squanti, Talon, bamboo, flair considered independent read as full length. Do the community has any suggestions?
Thanks

I haven't used these long-read RNA-seq tools so I could be wrong here. Transcript reconstruction combines reads from the same transcript and should be robust to partial transcripts. In addition, Nanopore direct RNA-seq produces a lot of partial transcripts. I tend to believe many transcript reconstruction tools can deal with them.
As to what tool to use, StringTie is a popular one and has been shown to be better than Flair. IsoQuant seems promising. Its developers argued it is more accurate than StringTie, Bambu, Talon and Flair. Squanti only works for iso-seq.

You could have a go using LAST. I've been developing and using a protocol for mapping cDNA reads that is based around LAST-mapping since 2017. The most recent RY model addition in LAST seems to clean up the majority of the random gene counting noise issues that we were having.
See the gene counting part of the protocol here. If you're working from direct/native RNA rather than cDNA (or don't care about strand-specific mapping), then the strand correction in previous steps is not needed:
https://dx.doi.org/10.17504/protocols.io.5qpvonn2bl4o/v17


text:
motif searching in a multi-fasta file
I am having trouble trying to search for motifs in a multi fasta file. I have used two techniques one gives me the name of the sequence where the motif is found but doesn't give me the the motif and its position. The other method does not return anything. Please see the code below, any assistance will be highly appreciated.
This is the protein sequence:
(Please note I have only put in the first 3 sequences, the actual file I am working with has 100 sequences)
>P1
MSDNGPQNQRNAPRITFGGPSDSTGSNQNGERSGARSKQRRPQGLPNNTASWFTALTQHG
KEDLKFPRGQGVPINTNSSPDDQIGYYRRATRRIRGGDGKMKDLSPRWYFYYLGTGPEAG
LPYGANKDGIIWVATEGALNTPKDHIGTRNPANNAAIVLQLPQGTTLPKGFYAEGSRGGS
QASSRSSSRSRNSSRNSTPGSSRGISPARMAGNGGDAALALLLLDRLNQLESKMSGKGQQ
QQGQXXTXKSAAEASKKPRQKRTATKAYNVTQAFGRRGPEQTQGNFGDQELIRQGTDYKH
WPQIAQFAPSASAFFGMSRIGMEVTPSGTWLTYTGAIKLDDKDPNFKDQVILLNKHIDAY
KTFPPTEPKKDKKKKADETQALPQRQKKQQTVTLLPAADLDDFSKQLQQSMSSADSTQA*

>P2
MSDNGPQNQRNAPRITFGGPSDSTGSNQNGERSGARSKQRRPQGLPNNTASWFTALTQHG
KEDLKFPRGQGVPINTNSSPDDQIGYYRRATRRIRGGDGKMKDLSPRWYFYYLGTGPEAG
LPYGANKDGIIWVATEGALNTPKDHIGTRNPANNAAIVLQLPQGTTLPKGFYAEGSRGGS
QASSRSSSRSRNSSRNSTPGSSRGISPARMAGNGGDAALALLLLDRLNQLESKMSGKGQQ
QQGQTVTKKSAAEASKKPRQKRTATKAYNVTQAFGRRGPEQTQGNFGDQELIRQGTDYKH
WPQIAQFAPSASAFFGMSRIGMEVTPSGTWLTYTGAIKLDDKDPNFKDQVILLNKHIDAY
KTFPPTEPKKDKKKKADETQALPQRQKKQQTVTLLPAADLDDFSKQLQQSMSSADSTQA*

>P3
MSDNGPQNQRNAPRITFGGPSDSTGSNQNGERSGARSKQRRPQGLPNNTASWFTALTQHG
KEDLKFPRGQGVPINTNSSPDDQIGYYRRATRRIRGGDGKMKDLSPRWYFYYLGTGPEAG
LPYGANKDGIIWVATEGALNTPKDHIGTRNPANNAAIVLQLPQGTTLPKGFYAEGSRGGS
QASSRSSSRSRNSSRNSTPGSSRGISPARMAGNGGDAALALLLLDRLNQLESKMSGKGQQ
QQGQTVTKKSAAEASKKPRQKRTATKAYNVTQAFGRRGPEQTQGNFGDQELIRQGTDYKH
WPQIAQFAPSASAFFGMSRIGMEVTPSGTWLTYTGAIKLDDKDPNFKDQVILLNKHIDAY
KTFPPTEPKKDKKKKADETQALPQRQKKQQTVTLLPAADLDDFSKQLQQSMSSADSTQA*

The first one is:
infile=open("sequence.fasta",'r')
out=open("Result.csv",'w')
pattern=re.compile(r"[YFWLIMVA]..[LMALVN].[AGSTCD]..[LAIVNFYMW]")
for line in infile:
  line = line.strip("\n")
  if line.startswith('>'):
     name=line
else:
   s=re.finditer(pattern,line)
   print('%s:%s' %(name,s))
   out.write('%s:\t%s\n' %(name,s))

This one above returns, the output below. It gives me the name of the sequence which is P100, but it does not give me (1) The motif found (2) the position of the motif. And the Results.csv file that it generates is blank.
#output
>P100:<callable_iterator object at 0x000001D0611B7AF0

I then tried a different technique:
infile=open("sequence.fasta",'r')
open=open("result.csv",'w')
pattern=re.compile(r"(P[A-Z]{2}P")
for line in infile:
  line=line.strip("\n")
  if line.startswith('>'):
     name=line
else:
   s = re.finditer(pattern,line)
   for match_obj in s:
       print(match_obj)
       print(match_obj.group())
  out.write('%s:\t%s\t%s\n' %(name, match_obj,match_obj.group()))

The one above does not return anything at all. the Results.csv file is also blank.
I'm still new in this python, and some of the techniques I used here I found on bioinformatics stackexchange.
The output I desire, is where I can get (1) Name of sequence, (2) Motif found and (3) position of motif (where it starts and ends) as shown in the example below in the context of every sequence in the multi fasta file
>P100
motif:
position:

Any format is fine as long as  I can get something similar to the format above. Could I also be assistet with   the indention as that cofuses me.

@terdon I have made the amendments and put in more detail

Once again: please don't make us debug your example code. You need to make sure that the code you give us actually executes so we don't waste time on syntax errors. The code you have given does not run. for one thing, it is missing the import re. For another, you are still using open as a variable name and then trying to write to out. All of these would give an error if you tried to execute the code. So please make sure you give us the actual, real code that you are executing!

Cross-post https://www.biostars.org/p/9545061/#9545129

can we use it with a file in stockholm format?

First of all, you can't do this this way. Fasta sequences can, and indeed yours do, span multiple lines so in order to find a pattern, you need to join those lines into one to make sure you don't miss any that are found across more than one line.
Next, as a general rule, it is a bad idea to hardcode output file names. That means the user cannot control the output, and that any existing files will be overwritten, so just print to standard output and redirect to a file:
foo.py > out_file

You also don't want to hardcode the name of the input file since that just makes the script harder to use. Instead, read the input file name as an argument. With all this in mind, here is a working version of your script:
#!/usr/bin/env python3
import re
import sys

def print_matches(seq_name, matches):
     if matches:
          for match in matches:
               print(name,match.start(), match.end(),match.group(),sep='\t' )

pattern = re.compile(r"[YFWLIMVA]..[LMALVN].[AGSTCD]..[LAIVNFYMW]")
seq = ""

try:
     with open(sys.argv[1],'r') as infile:
          for line in infile:
               line = line.strip("\n")
               if line.startswith('>'):
                    name=line.lstrip('>')
                    if seq != "":
                         matches = re.finditer(pattern, seq)
                         print_matches(name, matches)
                         seq = ""
               else:
                    seq += line
          else:
               matches = re.finditer(pattern, seq)
               print_matches(name, matches)

except IOError as e:
     print("Failed to open file %s with error %s" % (sys.argv[1], e))
     sys.exit(1)

Running that on your example gives:
$ foo.py sequence.fasta
P2  44  53  LPNNTASWF
P2  122 131 YGANKDGII
P2  209 218 MAGNGGDAA
P2  219 228 ALLLLDRLN
P2  304 313 AQFAPSASA
P2  391 400 VTLLPAADL
P2  410 419 MSSADSTQA
P3  44  53  LPNNTASWF
P3  122 131 YGANKDGII
P3  209 218 MAGNGGDAA
P3  219 228 ALLLLDRLN
P3  304 313 AQFAPSASA
P3  391 400 VTLLPAADL
P3  410 419 MSSADSTQA
P3  44  53  LPNNTASWF
P3  122 131 YGANKDGII
P3  209 218 MAGNGGDAA
P3  219 228 ALLLLDRLN
P3  304 313 AQFAPSASA
P3  391 400 VTLLPAADL
P3  410 419 MSSADSTQA

The output is 4 tab-separated columns per match with the sequence name, the match's start position, its end position and the actual sequence matched.

Read the sequences with Biopython, reverse translate the sequences and the pattern peptides, construct pattern regexes, search (find_matches()) for the patterns for each sequence in a loop using DNA Chisel.
Disclaimer: I'm the maintainer of DNA Chisel


text:
FlowSOM multi-step clustering
There are a few different commonly used clustering algorithms within the single-cell space, although Leiden seems to be the top choice these days. FlowSOM is a classic package for analyzing flow cytometry data. It has a two-step approach for clustering. First, it builds a self-organizing map (SOM) where cells are assigned to 100 grid points. Then, it uses ConsensusClusterPlus to meta-cluster the SOM codes. I have not seen other packages that take this approach. There is no publication for FlowSOM and the documentation does not explain this part. Why this two-step approach and why isn't this more common?

The two step approach will be used in any cluster algorithm (except see "network theory" below), you just may not see it. Clustering methods are very common in genomics, they are popular because they are very visual and will stack against a heatmap nicely. They have been replaced by dimensionality reduction such as PCA type analysis. For data mining PCA is a much more robust method for a number of reasons for example they overcome the problems listed as "The Bad" and "The Ugly" below. PCA type analysis is difficult to visualise, particularly when it steps into multi-dimensions and thats why developers like clustering ... easy for the developer to understand and easy for the use to visually appreciate, particularly against a heat map.
Any classical clustering The first step is to build a matrix and the second is to reorganise the relationships within the matrix using a pre-defined algorithm, unweighted pairwise group method with arithmetic mean (UPGMA).

The matrix in this case will be the grid.
The second step is the clustering method. The method used ConsensusClusterPlus by Wilkerson and Hayes (2010) is, 0

The consensus clustering (CC) method provides quantitative and visual stability evidence for estimating the number of unsupervised classes in a dataset.

All that means is it stacks a heat map onto a clustering tree. Here's the output of the paper and the 'workhorse' is the tree presented in "A",

Like all great movies, it has 3 faces,
The Good using k-means clustering and that is well respected method in supervised-learning and suspect it performs pretty well in unsupervised learning too. It's a good method.
The bad There could easily be artefacts in the deeper relationships of the tree and that would need a full empirical investigation. I know that because genetic diversity is often measured by trees and loads of attention has gone into modelling the process and CCP is assuming homogeneity in the rate process - which rarely hplds in biology. Thus, there will be situations where that falls over, however, what I strongly suspect is that the k-means methods is attempting to circumvent this situation, without doing so explicitly (if that makes any sense). I would need to study the algorithm to assess this and I would not use clustering in data mining/unsupervised learning.
The Ugly The basic problem with clustering and the reason it isn't used in genetic diversity *, is because there's loads of clustering algorithms, they give different answers and there's no way to distinguish which is which. Having said all that k-means clustering is very respectable.
The reason it's not common is that there has never been a consensus on the cluster algorithm and modellers prefer to avoid them because they are only used for data exploration/mining and not useful for hypothesis testing. However I do stress that k-means supervised learning is a very good method, but that is not the same as clustering (which is unsupervised learning).
* or at least if it is used its just a method to visual a matrix and a like Warning clustering! will appear somewhere in the legend.

Parameters and conclusion There is a parameter issue with the Figure above is tracking manually. It demonstrates k has a HUGE influence on the result. In real world supervised learning parameterisation is automated, never ever would it be done manually. The technique is called parameterisation and it's kinda important (so like the user can't use a parameter that gives an output that suits them!! [I've not said that BTW]). So you'd ask an algorithm above the k=x algorithm to find the optimal value of k.
If FlowSOM automates the value for k and then tells you ... thats good and I'd be happy with that. If it doesn't, you have to manually feed loads of values of k and provide an independent assessment on which was the best. That is not cool and you will have to perform a laborious manual calculation.

Network theory Leiden looks to be part of network theory e.g. Trag and van Eck (2019). That is a significant advantage over classical clustering, because network theory is being dumped over the top, i.e. the third step, which can yield real power over the data set. Personally I follow graph theory, but both have major potential and are heavily under-utilised. Its not really fair to compare a classical method with the power of network theory. Certainly when entering graph theory what you can leverage from the graph (network) is enormous. I personally, thought the whole area was being ignored (there's lots of applications), but maybe its time will yet come into the mainstream.

To @Pallie's comments of
dimensionality and distance. I do not disagree at all and is the reason why dimensionality reduction approaches such as PCA or multi-dimensional scaling are preferred.
RAM overhead Of course the calculation will be hard on RAM: if k was fixed at 2 or 3, it will run reasonably quickly, but might not be meaning. Any k parameterisation will be RAM intensive, but IMO thats the nature of big data calculation and cloud computing: you just pick a computer with loads of RAM or just fit more RAM. There will also be a risk of course of RAM bottlenecks on standard machines (where it just freezes). The key thing is whether the approaches are giving biologically significant different results and personally I would consider it empirically. These are different calculations and this is complex data.

There is a lot of great info in your response. We might be counting the steps differently. Your step 1 ("the matrix in this case will be the grid") I would consider to be step 2. Step 1 would be generating the grid. Why generate the SOM grid if you can just run ConsensusClusterPlus directly on the samples?

People have started using Leiden (a continuation of Louvain, optimizing graph partitions for modularity score) because A) Curse of dimensionality makes distances meaningless in 20k dimensions B) K-means algorithms take far  too long to run and far too much memory on single cell sized datasets (100k * 20k).

@burger yes we are counting the steps differently. I would count your second step as two steps within clustering. I would see "your step 1" as preparing the data set for step 2 (which is my original step 1). There could be transformations or standardisations in your step 1.  The later being very likely, almost certain and without it the results could subject to significant artefacts.

@Pallie, thanks for your reply, I've addressed both points above. I certainly agree with the first point.

I just wanted to add a bit regarding dimensionality. My example was using FlowSOM, which is designed for flow data, which was limited to less than 10 markers when it was first released. Thus, unlike scRNA-seq with thousands of genes, there is not much additional dimensionality reduction possible. That may explain some of the choices they made.

Thanks @burger understand, yes agreed.


text:
Identify SNP from genetic position
I am looking for a way to obtain SNP names from their coordinates (Chromosome:BP) - very much the reverse of the question answered here:
(https://stackoverflow.com/questions/20251612/map-snp-ids-to-genome-coordinates)
Example: I have identified that a SNP at position 64377836 on Chromosome 17 is associated with a certain factor in my analysis and would like to find out whether it has a SNP ID and what this ID is.
Manually, I can search databases for Chr17:64377836 and find rs668 - but I would like to find a way to do this automatically, preferably in R.
Is this possible in R - and if so, how?

in what format you already have your input data? Is it a table, VCF? plink?

In a table - but I could easily reformat it.

You can do this in R using the library biomaRt from Bioconductor. Here an example starting from a table with 3 columns: Chr, Start, End. Positions are assumed to be in GRCh38.
library(biomaRt)
snpMart = useEnsembl(biomart = "snps",
             dataset = "hsapiens_snp")

# assuming your table is tab-separated, with column names
myData <- read.table("myTable.txt", sep="\t", header=True)

coords <- apply(myData, 1, paste, collapse = ":")

# submit the query
getBM(attributes = c('refsnp_id', 'chr_name', 'chrom_start', 'chrom_end', 'allele'),
  filters = c('chromosomal_region'),
  values = coords,
  mart = snpMart)

Source


text:
In VEP annotation, how is the codon field interpreted?
After annotating with VEP a VCF file, we obtain different fields. One of them is called Codons which represents the affected codon in the transcript of the gene. Below is a screenshot of Insertions from a sample:

Variants and INDELs affect genes differently depending on the transcript we are looking at. My question is, for some notations for Insertions like in the screenshot, how can we interpret these results? Let's take the example of -/TC. Is TC inserted at the end of the previous codon? or the next codon of the transcript?
Below are the different representations I have in mind where | is the position of the indel and ^ is the possible insertion scenario of TC.
                     Before insertion               after insertion
                           v|
first scenario          xxx xxx xxx                   xxx TCx xxx

second scenario         xxx xxx xxx                   xxx xxx TCx
                            |  ^

My gut is telling me that this depends more on the strand I am looking at. If the strand is +, it is likely to be the second scenario. Otherwise, it would be the first scenario. Am I thinking correctly? or does this notation have another meaning? Or maybe VEP is not able its position?
Thanks for the help.


text:
Manually order the set fails in UpsetR
This question was also asked on Biostars
I have tried to use UpsetR to visualise the input file, which can be found here. How is it possible to make UpsetR accept an order defined in sets=c("Nlab", "NQLD", "Ngla", "Nsyl", "Ntom", "Ntab", "Natt", "Cann", "Stub", "Slyc", "Vvin", "Atha")?
with keep.order = T

with keep.order = F

> library("UpSetR")
> orthogroups_df<- read.table("orthogroups.GeneCount.tsv",  header=T, stringsAsFactors = F)
> #All species
> selected_species <- colnames(orthogroups_df)[2:(ncol(orthogroups_df) -1)]
> selected_species
 [1] "Atha" "Cann" "NQLD" "Natt" "Ngla" "Nlab" "Nsyl" "Ntab" "Ntom" "Slyc" "Stub" "Vvin"
> head(orthogroups_df)
  Orthogroup Atha Cann NQLD Natt Ngla Nlab Nsyl Ntab Ntom Slyc Stub Vvin Total
1  OG0000000    0    0  965    0    0    3    0    0    0    0    0    0   968
2  OG0000001    0    1    3    0    0  448    0    0    0    0    0    0   452
3  OG0000002    0    1  313    0    0  120    1    0    1    0    0    0   436
4  OG0000003    0   93   15   21   46   16   33   63   36   25   39   26   413
5  OG0000004    1   42    2   34  109    6    8  154   11    9    4    0   380
6  OG0000005    0    2   61    1   34   44   91   70   43   20    1    0   367
> ncol(orthogroups_df)
[1] 14
> orthogroups_df[orthogroups_df > 0] <- 1
> # we only show intersections of interest ,
> intersections=list(list(selected_species),
+                    list("NQLD", "Ngla", "Natt", "Nlab", "Nsyl", "Ntab", "Ntom"),
+                    list("Stub", "Slyc"),
+                    list("Atha", "Vvin"),
+                    list("Ntab", "Nsyl", "Ntom"),
+                    list("Nlab", "NQLD", "Ngla"),
+                    list("Nlab", "NQLD", "Nsyl"),
+                    list("Nlab", "Ngla", "Nsyl"),
+                    list("NQLD", "Nsyl", "Ngla"))
> upset(orthogroups_df,
+       text.scale = c(1.4),
+       sets=c("Nlab", "NQLD", "Ngla", "Nsyl", "Ntom", "Ntab", "Natt", "Cann", "Stub", "Slyc", "Vvin", "Atha"),
+       order.by = "freq",
+       keep.order=T,
+       intersections = intersections,
+       sets.x.label="Total number of orthogroups",
+       mainbar.y.label = "Number of orthogroups")

It's possible that upset is using factor levels for ordering, which is the default approach for ggplot2. If this is the case, then it can be fixed by converting the variable to a factor, and using an explicit level definition:
sets <- factor(sets, levels=sets)

For your specific example:
upSets <- c("Nlab", "NQLD", "Ngla", "Nsyl", "Ntom", "Ntab", "Natt", "Cann", "Stub", "Slyc", "Vvin", "Atha")
upSets <- factor(upSets, levels=upSets)
upset(orthogroups_df,
   text.scale = c(1.4),
   sets=upSets,
   order.by = "freq",
   keep.order=T,
   intersections = intersections,
   sets.x.label="Total number of orthogroups",
   mainbar.y.label = "Number of orthogroups")


text:
Input for --score in PLINK 1.9 when generationg PRS using GWAS summary statistics
I'm conducting a simulation and I need to obtain polygenic risk scores (PRS) using genome wide association studies (GWAS) summary statistics.
It is known that the GWAS will give BETA and OR for quantitative and binary phenotypes respectively.
When conducting the simulation for quantitative traits I used the following code to obtain PRS from GWAS performed:
./plink1.9 --bfile binary_file1_name --linear interaction --pheno pheno_file_name --covar cov_file_name --parameters 1 --allow-no-sex
./plink1.9 --bfile binary_file2_name --score plink.assoc.linear 2 4 7 header --out prs_quantitative

This is fine, where plink.assoc.linear file contains BETA.
When conducting simulation for binary traits I used the following code to obtain PRS from GWAS performed:
./plink1.9 --bfile binary_file1_name --logistic interaction --pheno pheno_file_name --covar cov_file_name --parameters 1 --allow-no-sex
./plink1.9 --bfile binary_file2_name --score plink.assoc.logistic 2 4 7 header --out prs_binary

I'm not sure whether we need to transform OR to BETA before inputting the file to --score in performing this for binary traits or, simply use the original plink.assoc.logistic which has OR (instead of BETA) to go ahead.
Simply, do we need to input file containing BETA values "only" for --score or either a file containing OR values or BETA values is fine?


text:
RNA folding at specific temperature with ViennaRNA in python
I am trying to get dot-bracket-notations of single stranded RNA via the ViennaRNA python package (https://pypi.org/project/ViennaRNA/) at different temperatures. I have read in the docs (https://www.tbi.univie.ac.at/RNA/ViennaRNA/doc/html/examples_c.html) that the folding prediction at different temperatures is possible with C..but nowhere it is described how to specify the temperature when using python. Does anyone know how to specify the temperature at which the RNA folding is predicted with the ViennaRNA package in python?
Thanks for the help!

I think that the "Global Variables" section of the docs describes this:

For the Python interface(s) SWIG places global variables of the C-library into an additional namespace cvar. For instance, changing the global temperature variable thus becomes

RNA.cvar.temperature = 25

Globals are somewhat surprising to me as a way of handling this. It seems like you should be able to specify this differently. Maybe I'm missing something?

Very cool, thanks.