text:
Approach for alignment-free phylogenetic tree generation using kmers
I have a general question about building phylogenetic trees without alignment:
I want to build a tree using an alignment-free approach by counting kmers. I have found a lot of information in the literature about this type of practice but I'm having trouble finding approaches that use Illumina short read sequences (which is what I have).
In my first attempt, I was able to generate a dendrogram by counting kmers and clustering them, but it did not group the sequences I had in the “correct/expected” way (these are known strains that are part of specific subclasses related to S. cerevisiae - in theory, this kmer clustering approach should’ve been able to group the strains based on their subclass pretty easily).
There are 21 samples in total with 5 specific subclasses. Does anyone know a way to either cluster or simply calculate the distance between kmers to generate a tree of 20+ Illumina short-read samples?
EDIT:
The main reason for doing this in this way is to eliminate any bias from the reference sequence and maintain all the alleles (most of the samples are tetraploid) so I really want to avoid collapsing the data with any type of preprocessing alignment or assembly.
With that being said I would want to use the whole-genome short-read sequences of these 21 samples to count the kmers across the genome for phylogenetic analysis. If it gets too complicated, I would be okay with assembling the reads first before kmer counting.
I used Kitsune to determine the optimal kmer size for the samples which comes out to 15mer. I was able to count 15mers with all the samples using three methods: pure raw reads, assembling the reads with SPAdes, and downsampling the raw reads to 5x coverage. Using Jellyfish and custom scripts I wrote, which didn't end up being too taxing on my computer/server.
Currently, I have 21 '.tsv' files that list the kmer and the kmer frequency for each individual sample and I have one giant matrix that combines that data together so I have the kmer and the kmer count for each sample in one file:

sample_1
sample_2
sample_3
...
sample_21

AAAAA
123432
345432
56765
...
55677

AAAAG
1023
0
350655
...
1

AAAAT
462
741
2136
...
121

AAAAC
342
123
3453
...
12345

------------
----------
----------
---------------

I have looked briefly at kWIP. Ideally, I would use their distance calculation function to generate a .dist file for my giant matrix - I'll look into what inputs can be used for that function.

I added more information about the specifics of the file I'm working and the reasoning if that helps!

A very simple solution is to use a correlation r or r2. This establishes a pairwise matrix which is then resolved under neighbour-joining (don't use UPGMA or similar). MegaX will import a pairwise matrix and produce a tree in neighbour-joining. You could even produce a bootstrapped solution (that would require custom code) - bit tricky.
I definitely recommend neighbour-joining for the tree
MegaX is a good implementation of neighbour-joining (it was originally Nei's idea and that was his package). The reason for this are complex and happy to address this in a separate question.
Algorithm Lets say r or r2 with Sample 1 versus Sample 2 ... the correlation value goes in one half of a triangular matrix with zeros running down the centre of the matrix.  You then proceed for every pairwise comparison until you get a triangle with Sample 1, Sample 2, Sample 3, etc ... along the top and  Sample 1, Sample 2, Sample 3, etc  ... along the side. That will import to MegaX (its a long while since I've used it BTW).
Distance Methods If I was performing this I wouldn't use  r or r2: I would use either Jaccard, Euclidean or Cosine (very trendy right now) distances. Cosine is trendy because taps into tensors (Tensorflow and all that). Scipy will almost certainly do the distance methods listed here for your data and output a nice triangular matrix.

Personally if I was coding this I would absolutely perform a bootstrap, but this is a lot of coding.
Issues The issue with distance metrics without bootstrapping is that they give different solutions and there's no way identify which tree is the best tree. Bootstrapping helps a lot here and values > 80% are considered robust topologies. I've never personally known two distance metrics produce an incongruence for the same data under bootstrapping. The distance matrix issue was why everyone moved (and rigidly stuck with) formal models via maximum likelihood or Bayes - but these require alignments.
Summary Think you can go alignment free? Well yeah you can, but remember there's no such thing as a free lunch ;-)

Okay, I love this thank you so much! My main takeaway is to get the distance matrix from the giant kmer count matrix I currently have and then do some neighbor-joining to get a tree which I can visualize with another tool later.

Yep - providing the tree method is called  "neighbor-joining" (not UPGMA anything like that) ... again if you want that explaining - long explanation - please submit another question.


text:
How to calculate average from a column when consecutive cells are similar in different columns?
I wonder how can I calculate the average of column 7 (WPS_win) when consecutive cells in column 2, column 3, and column 1 are the same. The file is in bed format (tab-separated format). The original file is very large. The first few lines of the file look like this,
chrom   start_pr    end_pr  width_pr    start_win   end_win WPS_win .overlap
chr1    981918  984117  2200    981900  982000  -5  82
chr1    981918  984117  2200    982000  982100  -1  100
chr1    981918  984117  2200    982200  982300  -2  100
chr1    981918  984117  2200    982300  982400  2   100
chr1    981918  984117  2200    982400  982500  -3  100
chr1    981918  984117  2200    982500  982600  -1  100
chr1    981918  984117  2200    982600  982700  -3  100
chr1    981918  984117  2200    982700  982800  1   100
chr1    981918  984117  2200    982800  982900  -5  100
chr1    981918  984117  2200    982900  983000  0   100
chr1    981918  984117  2200    983000  983100  0   100
chr1    981918  984117  2200    983100  983200  -5  100
chr1    981918  984117  2200    983200  983300  -1  100
chr1    981918  984117  2200    983300  983400  -3  100
chr1    981918  984117  2200    983400  983500  1   100
chr1    981918  984117  2200    983500  983600  -2  100
chr1    981918  984117  2200    983600  983700  1   100
chr1    981918  984117  2200    983700  983800  -3  100
chr1    981918  984117  2200    983800  983900  1   100
chr1    981918  984117  2200    983900  984000  -4  100
chr1    981918  984117  2200    984000  984100  -1  100
chr1    981918  984117  2200    984100  984200  -2  17
chr1    999138  1001337 2200    999100  999200  1   62
chr1    999138  1001337 2200    999200  999300  -1  100
chr1    999138  1001337 2200    999700  999800  -1  100
chr1    999138  1001337 2200    999800  999900  -2  100
chr1    999138  1001337 2200    999900  1000000 -2  100
chr1    999138  1001337 2200    1000000 1000100 1   100
chr1    999138  1001337 2200    1000100 1000200 -2  100
chr1    999138  1001337 2200    1000200 1000300 0   100
chr1    999138  1001337 2200    1000300 1000400 -1  100
chr1    999138  1001337 2200    1000400 1000500 -2  100
chr1    999138  1001337 2200    1000500 1000600 -2  100
chr1    999138  1001337 2200    1000600 1000700 2   100
chr1    999138  1001337 2200    1000700 1000800 -1  100
chr1    999138  1001337 2200    1000800 1000900 -2  100
chr1    999973  1002172 2200    999900  1000000 -2  27
chr1    999973  1002172 2200    1000000 1000100 1   100
chr1    999973  1002172 2200    1000100 1000200 -2  100
chr1    999973  1002172 2200    1000200 1000300 0   100
chr1    999973  1002172 2200    1000300 1000400 -1  100
chr1    999973  1002172 2200    1000400 1000500 -2  100
chr1    999973  1002172 2200    1000500 1000600 -2  100
chr1    999973  1002172 2200    1000600 1000700 2   100
chr1    999973  1002172 2200    1000700 1000800 -1  100
chr1    999973  1002172 2200    1000800 1000900 -2  100
chr1    999973  1002172 2200    1001600 1001700 -1  100
chr1    999973  1002172 2200    1001700 1001800 -2  100
chr1    999973  1002172 2200    1001800 1001900 1   100
chr1    999973  1002172 2200    1001900 1002000 -1  100

In the above table, the first 22 rows of cols 1, 2, and 3 are the same. Therefore, The expected output would look like this,
chrom   start_pr    end_pr  Mean_WPS
chr1    981918  984117  -1.590909091
chr1    999138  1001337 -0.857142857
chr1    999973  1002172 -0.857142857

I'm open to both R and Bash. Any suggestion to overcome this issue would be heartily appreciated.

library(dplyr)
cols_of_interest <- c("chrom", "start_pr", "end_pr")

test_data |>
  group_by(across(all_of(cols_of_interest))) |>
  summarize(Mean_WPS = mean(WPS_win))

gives the following with your test data above:
# A tibble: 3 × 4
# Groups:   chrom, start_pr [3]
  chrom start_pr  end_pr Mean_WPS
  <chr>    <int>   <int>    <dbl>
1 chr1    981918  984117   -1.59
2 chr1    999138 1001337   -0.857
3 chr1    999973 1002172   -0.857

Ref for the code used above

I use bedtools merge with a little trick to combine the first three columns into a single one to "mock" each of the unique intervals into one "chromosome" so bedtools easily knows that these are the same:
cat example.bed \
| awk 'NR > 1 {OFS="\t"; print $1"--"$2"--"$3"--", 1, 1, $7}' \
| bedtools merge -i - -c 4 -o mean \
| awk 'BEGIN {OFS="\t"; print "chrom", "start_pr", "end_pr", "Mean_WPS"} {split($1,x,"--"); print x[1], x[2], x[3], $4}'

chrom   start_pr        end_pr  Mean_WPS
chr1    981918  984117  -1.590909091
chr1    999138  1001337 -0.8571428571
chr1    999973  1002172 -0.8571428571

First awk command: Skips first line, combines first three columns into one and prints this together with a single position for start and end plus the score.
bedtools command: Merges the identical intervals, taking the mean of the fourth column.
Second awk command: Splits the first column again into the former chrom-start-end and puts the average score in column 4, and also adds a header line.
Memory footprint of this should be minimal.

Side note: If one copies the example data from OP then simply do a tr -s " " | tr " " "\t" on it to remove the many whitespaces and replace by a single one, then replace whitespace by tab to create a proper tsv file for parsing.

source: Stack Exchange

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

sort: votes


text:
which SNP effect software, such as VEP, annovar or snpEff, allows the user to submit their own annotations in gtf format?
I have a use case where I want to annotate the variants from running Oxford Nanopore's Medaka on a very simple reference, a single sequence starting on an ATG and finishing on a stop codon (no introns).
Which of the variant effect prediction tools would be the easiest to use to programmatically create the .gtf annotation file and build a reference upon which one can then call the variant annotation, and get the different types of variant consequences on the vcf file (synonymous, non-synonymous with the pair of aa's involved, frame-shifts, etc.)
So far I have looked at:

snpEff: it seems to have a way of building the reference, although it's not straightforward to match the built java .jar with the one in the system.

annovar: I couldn't find the documentation on how to supply your own annotations, only how to download the human GRCh38 or equivalents.

VEP: I couldn't find the documentation on how to supply your own annotations, only how to download the human GRCh38 or equivalents.

Any ideas?

I managed to solve the issue by using a slightly different approach to @Ram RS, which is bcftools csq instead of bcftools annotate. The formatting is via GFF3, and the tool is picky about how to format it, but it turns out Torsten @tseemann had already gone through the pain of matching what the tool expects with the GFF3 flavour, and then it works:
https://github.com/samtools/bcftools/issues/530#issuecomment-377637623
An example gff3 for a single CDS of length=1463 is:

seq1    geneX   gene    1       1463    .       +       .       ID=gene:gene1;Name=gene1;biotype=protein_coding
seq1    geneX   mRNA    1       1463    .       +       .       ID=transcript:transcript1;Name=transcript1;Parent=gene:gene1;biotype=protein_coding
seq1    geneX   CDS     1       1463    .       +       0       ID=CDS:cds1;Name=CDS1;Parent=transcript:transcript1;biotype=protein_coding

Then call it with:
bcftools csq -p m -f $fastaref -g $gff3 $inputvcf -o $outvcf

Thank you all for commenting.

All of them have ways you can add your own annotation dataset but I've found the easiest way to be to use bcftools annotate. You can give it a tab-delimited (GTF is tab-delimited, but might need some pre-processing to get relevant values from column-9) or VCF file and have it pick specific columns/attributes, match the locus by any subset of CHR, POS, REF, ALT and add the annotation in a variety of custom manners. vcfanno is another option.

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


text:
Find protein in DNA sequences with arbitrary encoding and possible frameshifts
I have a dataset with lots of DNA sequences (~10 M seqs, each ~6 KBp long).
I want to detect presence of a fixed given protein (~600 Bp / 200 aa) within each sequence and, if detected, obtain the frame and position of the protein within the sequence. However, the protein may be coded arbitrarily and differently in each sequence. In addition, I need to account for frameshifts induced by sequencing errors.
I've tried to use tfasty36 from the FASTA36 package. This is basically a tblastn equivalent (protein query search against DNA database in all frames). Due to the statistical filtering built into these tools, I've come to the conclusion that tfasty36 and tblastn are not the right tools for this job.
tfasty36 didn't work because it requires fiddling with statistical tolerances to show all matches. Also, the matches are quite fragmented (I guess, it uses heuristics instead of dynamic programming). For tblastn, building the database seems to be unnecessary (given that it may have the same problems as tfasty36 afterwards).
I've looked into alignment tools, but couldn't find any that can align a protein sequence to a DNA sequence while taking frameshifts into account.
How can I solve this problem?

If you need to model introns and splicing, then you should use a tool that is designed for that. One of the limitations of blast and its cousins is that they only look at sequence similarity and do not attempt to join HSPs into possible transcript models.
I haven't worked on this in a long time, but back when I did, the programs I used were:

exonerate

Exonerate is a generic tool for pairwise sequence comparison. It allows you to align sequences using a many alignment models, either exhaustive dynamic programming or a variety of heuristics.

In this case, you would use your protein sequence as the query and the set of  DNA sequences as the target and try to align using the protein2genome model. Something like this:
  exonerate --model protein2genome protein.pep dna.fa > out

This will result in (many) alignments of your query protein agains the various genomic sequences. Here's an example I had handy:
$ exonerate --model protein2genome -n 1 frog.histone.pep c.elegans.chrI.fa
Command line: [exonerate --model protein2genome -n 1 frog.histone.pep c.elegans.chrI.fa]
Hostname: [oregano]

C4 Alignment:
------------
         Query: frog
        Target: CM041006.1 Caenorhabditis elegans strain DL226 chromosome I, whole genome shotgun sequence:[revcomp]
         Model: protein2genome:local
     Raw score: 732
   Query range: 32 -> 300
  Target range: 3039620 -> 3034073

      33 : LeuGlnLysLysLeuGluGluLeuGluLeuAspGluGlnGlnArgLysArgLeuGl :      51
           |||:!!:!!:!!|||:!!:!!|||! !|||   |||! !|||! !:!!||||||..
           LeuGluGluGlnLeuLysLysLeuGlyLeuThrGluProGlnThrGlnArgLeuSe
 3039620 : CTAGAAGAACAGCTGAAAAAGCTCGGGCTCACGGAGCCACAAACTCAACGGCTCAG : 3039566

      52 : uAlaPheLeuThrGlnLysGlnLysValGlyGluLeuLysAspAspAspPheGluL :      70
           .! !||||||  !   |||:!!   :!!  !||||||!..!!:|||   .!.:!!!
           rGluPheLeuGlnValLysGluGlyIleLysGluLeuSerGluAspMetLeuGlnT
 3039565 : CGAATTCTTGCAGGTCAAGGAGGGCATCAAGGAGTTGAGTGAGGATATGCTACAGA : 3039509

      71 : ysValSerGluLeuGlyAlaGlyAsnGlyGlyValValPheLys  >>>> Targe :      85
             ! !.!!|||||||||  !||||||||||||||||||  !|||
           hrGluGlyGluLeuGlyHisGlyAsnGlyGlyValValAsnLys++
 3039508 : CTGAAGGAGAACTTGGCCATGGAAATGGAGGAGTTGTTAATAAGgt.......... : 3039462

      86 : t Intron 1 >>>>  ValSerHisLysProThrSerLeuIleMetAlaArgLys :      97
           55 bp              !  !|||::!  !|||.!!:!!|||||||||||||||
                          ++CysValHisArgLysThrGlyValIleMetAlaArgLys
 3039461 : ...............agTGTGTTCATCGAAAAACCGGTGTAATTATGGCTCGGAAG : 3039373

      98 : LeuIleHisLeuGluIleLysProAlaIleArgAsnGlnIleIleArgGluLeuGl :     116
           |||:!!||||||||||||||||||:!!:!!|||.!.||||||:!!!:!||||||
           LeuValHisLeuGluIleLysProSerValArgGlnGlnIleValLysGluLeuAl
 3039372 : CTTGTACACCTCGAGATCAAGCCGTCAGTCCGACAACAGATCGTCAAAGAGCTCGC : 3039316

     117 : nValLeuHisGluCysAsnSerProTyrIleValGlyPheTyrGlyAlaPheTyrS :     135
            |||||||||:!!||||||||||||!:!||||||||||||||||||||||||   .
           aValLeuHisLysCysAsnSerProPheIleValGlyPheTyrGlyAlaPheValA
 3039315 : CGTACTTCACAAGTGCAATAGTCCATTCATTGTTGGATTCTACGGGGCTTTTGTGG : 3039259

     136 : erAspGlyGluIleSerIleCysMetGluHisMetAspGlyGlySerLeuAspGln :     153
           .!:!!..!!!:||||||||||||||||||:!!|||||||||  !|||||||||  !
           spAsnAsnAspIleSerIleCysMetGluTyrMetAspGlyLeuSerLeuAspIle
 3039258 : ACAATAACGATATATCAATTTGTATGGAGTATATGGATGGGCTTTCACTTGATATA : 3039205

     154 : ValLeuLysLys  >>>> Target Intron 2 >>>>  AlaGlyLysIlePro :     162
           ||||||||||||           3837 bp           !.!|||::::!!|||
           ValLeuLysLys++                         ++ValGlyArgLeuPro
 3039204 : GTTTTGAAAAAGgt.........................agGTTGGCCGCCTCCCG : 3035341

     163 : GluLysIleLeuGlyLysValSerIleAlaValIleLysGlyLeuThrTyrLeuAr :     181
           ||||||.!!:!!|||::!:!!|||:!!||||||:!!:::|||||||||||||||!:
           GluLysPheValGlyArgIleSerValAlaValValArgGlyLeuThrTyrLeuLy
 3035340 : GAAAAGTTTGTCGGACGAATCTCGGTTGCAGTTGTCCGTGGATTGACGTACTTGAA : 3035284

     182 : gGluLysHisLysIleMetHisArg{A}  >>>> Target Intron 3 >>>>  :     190
           !!!::!!  !||||||:!!||||||{|}            50 bp
           sAspGluIleLysIleLeuHisArg{A}++                         +
 3035283 : AGATGAAATTAAGATTTTGCATAGA{G}gt.........................a : 3035208

     191 :  {sp}ValLysProSerAsnIleLeuValAsnSerArgGlyGluIle--LysLeuC :     207
            {||}|||||||||||||||!!:||||||||||||!.!|||||||||##   ||||
           +{sp}ValLysProSerAsnMetLeuValAsnSerAsnGlyGluIle##---LeuC
 3035207 : g{AT}GTTAAACCATCGAATATGCTTGTCAATAGTAACGGAGAAATTAA---TTAT : 3035157

     208 : ysAspPheGlyValSerGlyGlnLeuIleAspSerMetAlaAsnSerPheValGly :     225
           ||||||||||||||||||||..!|||||||||||||||||||||||||||||||||
           ysAspPheGlyValSerGlyMetLeuIleAspSerMetAlaAsnSerPheValGly
 3035156 : GCGATTTTGGAGTCTCTGGAATGTTGATTGATTCGATGGCCAACTCATTTGTCGGA : 3035103

     226 : ThrArgSerTyrMetSer  >>>> Target Intron 4 >>>>  ProGluArg :     234
           |||||||||||||||:!!            799 bp           |||||||||
           ThrArgSerTyrMetAla++                         ++ProGluArg
 3035102 : ACGAGGAGTTATATGGCGgt.........................agCCCGAACGA : 3034277

     235 : LeuGlnGlyThrHisTyrSerValGlnSerAspIleTrpSerMetGlyLeuSerLe :     253
           |||  !|||:!!||||||!!!:!!..!|||||||||||||||.!.|||||||||||
           LeuThrGlySerHisTyrThrIleSerSerAspIleTrpSerPheGlyLeuSerLe
 3034276 : CTCACAGGATCCCACTATACAATTTCATCGGATATCTGGTCATTTGGATTATCACT : 3034220

     254 : uValGluMetAlaIleGlyArgTyrProIleProProProAspAlaLysGluLeuG :     272
           |||||||:!:  !|||||||||||||||:!:||| !!|||..!  !  !|||! !!
           uValGluLeuLeuIleGlyArgTyrProValProAlaProSerGlnAlaGluTyrA
 3034219 : TGTCGAACTTCTTATTGGACGATATCCAGTGCCTGCACCGTCTCAAGCAGAATATG : 3034163

     273 : luLeuIlePheGlyCysSerVal<->GluGlyAspProAlaSerSerGluLeuAla :     289
               !!!:|||..!   :!!! !   |||  !!!:! !|||..!|||     !! !
           laThrMetPheAsnValAlaGluAsnGluIleGluLeuAlaAspSerLeuGluGlu
 3034162 : CTACAATGTTTAATGTGGCTGAAAATGAGATTGAACTGGCTGATAGTCTTGAAGAG : 3034109

     290 : ProArgProArgProProGlyArgProIleSer :     300
           |||!.!  !!.!||||||.!!!.!|||   |||
           ProAsnTyrHisProProSerAsnProAlaSer
 3034108 : CCCAACTATCATCCGCCGAGCAATCCGGCGTCA : 3034074

vulgar: frog 32 300 . CM041006.1 3039620 3034073 - 732 M 52 156 5 0 2 I 0 51 3 0 2 M 73 219 5 0 2 I 0 3833 3 0 2 M 32 96 S 0 1 5 0 2 I 0 46 3 0 2 S 1 2 M 14 42 F 0 2 G 1 0 M 26 78 5 0 2 I 0 795 3 0 2 M 48 144 G 0 3 M 21 63
-- completed exonerate analysis

In addition, exonerate lets you control the output you want. Check out the --ryo (roll your own output) option in man exonerate.

genewise

GeneWise compares a protein sequence to a genomic DNA sequence, allowing for introns and frameshifting errors.

Unlike exonerate, genewise only does protein to DNA alignment, but like exonerate it models introns and splicing. 20 years ago, when I did this every day, my impression was that genewise was slower but slightly more sensitive than the default exonerate options (exonerate can be tweaked to be made more sensitive). So I would usually try exonerate first and fall back to genewise if it failed to find a good match.


text:
Error importing Maximum Credibility Tree from BEAST into R using treeio's read.beast()
I've encountered an issue that I'm hoping someone can help with. I constructed a phylogenetic tree using BEAST1 and then extracted the Maximum Credibility Tree using treeannotator. I'm trying to import it into R for editing, but I run into an error when using the read.beast() function from the treeio package:
MCC_tree <- read.beast("/home/MyPC/Downloads/MCC_beast.tre")
Error in FUN(X[[i]], ...) :
  numbers of left and right parentheses in Newick string not equal

Interestingly, when I attempt to visualize the tree using figtree, it works fine. Likewise, importing into R using read.nexus() works without issues. However, I want to use read.beast() specifically because I need to annotate the "posterior", which I can't do when I import it using read.nexus().
Has anyone encountered a similar issue or have any suggestions on how to resolve this?
Thank you in advance!

if you read it into FigTree then in the export in Nexus and perform a "save all trees"  it should then read into read.beast().
Use the beast file output before you performed tree annotator operation.

1 or 2. should be read directly into read.beast().
I would try option 2 first.


text:
A function to produce M fragments of DNA with length of I
I want a script in R, which creates M fragments of DNA with a length of I but with GC content of >70%
I know this gives something but not exactly what I am looking for
require(Biostrings)

set = NULL
for (i in 1:4)
  set = c(set, paste(sample(DNA_ALPHABET[1:4], 3, replace=T), collapse=""))
rep(set,2)

Any idea please?

Why R though? If the need is to simulate sequences, why do you need to do it R - why not look for the best way to do it using any language/tool?

Also by this point you should know better than to use set as a variable name.

There are different ways to approach the GC bias, the easiest is probably to decide beforehand how many "A/T" and how many "G/C" we need in each sequence, then choose them randomly.
# set global parameters
M <- 7
I <- 15

# for each sequence, choose a proportion of GC
prop_gc <- runif(M,
                 min = .7,
                 max = 1)

# randomly choose bases, ensure we have a number that gives us at least `prop_gc`
bases <- lapply(prop_gc,
                \(.prop_gc){
                  c(
                    sample(c("A","T"),
                           size = floor(I*(1 - .prop_gc)), replace = TRUE),
                    sample(c("G","C"),
                           size = ceiling(I*.prop_gc), replace = TRUE)
                  )
                })

# now we can simply take our bases and combine them (in a random order)
sequences <- sapply(bases,
                    \(.nt) sample(.nt) |> paste(collapse = ""))

sequences
#> [1] "CGGCTTCCGCGCGGGT" "CTACCCACGGCGCCGG" "GGGCGGTGTATGTGCC" "GGCCCGGGTGCCCAGC"
#> [5] "GCGGGCGCGGGCGGGT" "CACCCCCTCCCGGGAC" "GTCGAAGCGCGCCAGC"

Created on 2023-10-12 with reprex v2.0.2
Or same thing with a for loop (since you gave it as a tag):
M <- 7
I <- 15

sequences <- character(length = M)
for(i in seq_len(M)){
  prop_gc <- runif(1,
                   min = .7,
                   max = 1)

  bases <- c(
    sample(c("A","T"),
           size = floor(I * (1 - prop_gc)), replace = TRUE),
    sample(c("G","C"),
           size = ceiling(I * prop_gc), replace = TRUE)
  )

  sequences[[i]] <- sample(bases) |> paste(collapse = "")

}
sequences
#> [1] "CCCCGGTTGCGTTTCC" "CCTCGCTGCGCGTGGC" "CAGTCTCCCGCGGGGC" "GGGCCCGCCGAAAGGC"
#> [5] "GCGCGGGGCCGGTCAG" "CGGGCGCCCTGGACAC" "ACTCTGCGGGGGAGTC"

Created on 2023-10-12 with reprex v2.0.2
And we can doublecheck the results:
M <- 7000
I <- 500

[...]

# check
gc_prop <- function(x){
  toupper(x) |>
    strsplit("") |>
    sapply(\(.char) mean(.char == "C" | .char == "G"))
}

sapply(sequences, nchar) |> table()
#>
#>  500
#> 7000

min(gc_prop(sequences))
#> [1] 0.702

Created on 2023-10-12 with reprex v2.0.2


text:
MegaX Pairwise Matrix Error
I have a pairwise matrix that I've manually formatted into a ".mega" file following the guide in the MegaX documentation but when I go to import it into the application it gives this error:
Unable to find the command terminator '
' (in line 48)

I don't understand what the terminator is supposed to be I can't find anything in the documentation or online. This is what my file looks like:
#mega
TITLE: Lower-left triangular matrix

#sample1
#sample2
#sample3
#sample4
#sample5
#sample6
#sample7
#sample8
#sample9
#sample10
#sample11
#sample12
#sample13
#sample14
#sample15
#sample16
#sample17
#sample18
#sample19
#sample20
#sample21
#sample22

6779.024
4010.596    5085.012
4444.925    3965.719    3597.515
3824.67     4435.878    2140.942    2855.545
2698.035    5877.724    4274.36     3721.735    3562.716
5746.71     4991.678    3103.77     4587.394    3318.78     6141.058
2898.099    7345.178    5461.126    5000.322    4846.293    2354.302    7366.433
4667.862    2949.764    3986.194    2832.598    3118.544    3551.237    5007.972    5056.938
3956.69     5110.169    2489.003    3530.347    2310.696    4285.076    3286.188    5301.894    3994.313
4086.841    5636.588    2851.673    3681.433    3078.399    4197.933    4416.542    5139.341    4347.211    3250.879
3595.868    6237.556    3924.761    3907.091    3624.064    3081.475    5770.594    3985.626    4223.117    4013.377    3136.24
4574.567    5763.568    2420.15     4433.579    3203.997    5015.929    3445.061    6063.541    4820.971    3131.669    2940.224    4132.314
4621.723    3265.735    3518.476    2528.629    2685.442    3581.311    4361.296    5144.475    2322.969    3536.469    3942.852    4218.651    4431.758
4914.659    3750.498    2719.845    3117.886    2456.598    4411.935    3294.878    5908.595    3261.36     3157.398    3657.63     4540.83     3583.543    2199.326
5662.093    3090.623    3190.918    3289.562    2772.944    5317.999    2559.264    6786.94     3426.171    3247.437    4233.766    5294.141    3960.851    2656.578    2155.925
5559.707    3073.796    4997.709    3838.601    4294.455    4486.788    5728.029    5962.611    2303.294    5051.106    5270.304    5020.468    5639.198    3504.092    4111.158    4282.95
3906.494    4028.613    3618.28     2618.348    2835.36     2743.097    4941.507    4259.628    2191.873    3599.004    3847.302    3506.872    4432.612    2061.23     2947.972    3583.198    3342.373
6582.561    1663.359    4830.862    3730.956    4161.803    5594.268    4860.298    7160.793    2675.27     4859.183    5448.895    5895.575    5537.159    3011.78     3374.505    2881.87     2979.42     3464.931
4706.715    4004.427    4618.035    3537.473    3842.462    3773.964    5703.134    5002.92     2263.311    4592.613    4872.922    4301.443    5290.578    3530.202    4175.393    4494.308    2427.841    2677.65     3599.691
5703.386    2214.959    4515.044    3273.191    3766.439    4681.361    5022.943    6197.46     1837.686    4544.365    5026.108    5096.364    5257.581    2833.977    3419.825    3267.254    2034.999    2778.249    1661.474    2332.169
6689.288    11865.223   8990.382    9113.452    8897.297    7007.562    10983.956   6142.025    9502.64     9014.072    8056.924    6431.134    8821.925    9843.732    10166.579   11038.532   10009.974   8748.955    11573.491   8952.942    10488.639

EDIT:
I copied the contents of the .meg file over to an Excel sheet to make it easier to see but I changed the #Mega to #mega and I made sure to only include the correct elements (there are 22 samples in total, only 21 values present in the first line of the matrix). I also made sure to keep a line between the last sample name and the first value of the matrix. However, it's still giving me the same error.
Is there something I'm supposed to add at the end of the file to make it read the file correctly? I'm  trying to import this distance matrix so I can generate a NJ tree. I'm not sure what I'm doing differently than what is mentioned here.

Criteria, matrix. There must be a blank line between the last sample name and the start of the matrix. Also better to use #mega.
The tree from the matrix in the question is.

It's simply the carriage returns must by in Windows format \r (even if its Linux or Mac) and not Mac or Linux \n. In Mac the converter is BBEdit.

Okay I edited the original post again with what I have - hopefully this helps. I changed it to a lower-left matrix and I was able to make a plot using this matrix with R but I want to see how it works with MegaX

Okay it's done. Please do remember to upvote/accept this and the original answer to your first question from where this question originates.

Oh my goodness such a simple solution - thank you so much! I downloaded BBEdit and converted the file to a Windows (CRLF) .meg file and that worked!

Thanks @rimo good to know it worked. Mega was originally written as a Windows program - thats why I suspected the carriage return format. When a developer changes platform there's always "carry-over" from the previous code.


text:
Converting VCF format to text for use with PLINK and understanding column mapping
I successfully completed Nature PRS tutorial, which is based on PLINK. Turning to my real data, I downloaded ukb-d-20544_1.vcf.gz. Now I'm facing the problem that I seem to be unable to use it in PLINK or find the correct data format to download at all, and I am a bit confused about the correct general summary statistics-download procedure.
1. Convert .vcf format to a .txt file, PLINK can understand
In the tutorial they use a .txt.gz base data file aka summary statistics that looks like this:
CHR BP  SNP A1  A2  N   SE  P   OR  INFO    MAF
1   756604  rs3131962   A   G   388028  0.00301666  0.483171    0.997886915712657   0.890557941364774   0.369389592764921
1   768448  rs12562034  A   G   388028  0.00329472  0.834808    1.00068731609353    0.895893511351165   0.336845754096289
...

Having this .vcf.gz (please see excerpt at the very bottom), I unzipped and, in R, used MungeSumstats::format_sumstats("ukb-d-20544_1.vcf", ref_genome = "GRCh37")[1], to get something out of it similar to what they use in the .txt file in the tutorial. However, the column names won't match exactly. I know that I can use BETA instead of OR, and FRQ might correspond to MAF, but I can't seem to find the imputation information score INFO. Here is what MungeSumstats::format_sumstats produced:
SNP CHR BP  A1  A2  END FILTER  FRQ BETA    SE  LP  N   P
rs12238997  1   693731  A   G   693731  PASS    0.116734    0.00058531  0.000756397 0.357493    117706  0.439042942120051
rs371890604 1   707522  G   C   707522  PASS    0.0983981   0.000793177 0.000848559 0.456023    117706  0.349926634613021
...

I need to have the columns corresponding to MAF and INFO as per the tutorial.
The next step in the tutorial involves filtering based on specific criteria, such as MAF and INFO values. However, I can't proceed because these columns are missing.
I would like to primarily use the terminal and only turn to R for statistical analysis. Could you guide me on the correct approach to convert the .vcf file to .txt so that I can continue following the tutorial with PLINK?
2. General approach
In addition to addressing the conversion issue, I'm also seeking guidance on the appropriate data sources for genetic analysis. I'm curious about whether OpenGWAS is a reliable source for publication or if there are other recommended sources like dbGaP. Additionally, what is the ideal data format to obtain?

The beginning of the received .vcf file looks like this:
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
##FORMAT=<ID=ES,Number=A,Type=Float,Description="Effect size estimate relative to the alternative allele">
##FORMAT=<ID=SE,Number=A,Type=Float,Description="Standard error of effect size estimate">
##FORMAT=<ID=LP,Number=A,Type=Float,Description="-log10 p-value for effect estimate">
##FORMAT=<ID=AF,Number=A,Type=Float,Description="Alternate allele frequency in the association study">
##FORMAT=<ID=SS,Number=A,Type=Float,Description="Sample size used to estimate genetic effect">
##FORMAT=<ID=EZ,Number=A,Type=Float,Description="Z-score provided if it was used to derive the EFFECT and SE fields">
##FORMAT=<ID=SI,Number=A,Type=Float,Description="Accuracy score of summary data imputation">
##FORMAT=<ID=NC,Number=A,Type=Float,Description="Number of cases used to estimate genetic effect">
##FORMAT=<ID=ID,Number=1,Type=String,Description="Study variant identifier">
##META=<ID=TotalVariants,Number=1,Type=Integer,Description="Total number of variants in input">
##META=<ID=VariantsNotRead,Number=1,Type=Integer,Description="Number of variants that could not be read">
##META=<ID=HarmonisedVariants,Number=1,Type=Integer,Description="Total number of harmonised variants">
##META=<ID=VariantsNotHarmonised,Number=1,Type=Integer,Description="Total number of variants that could not be harmonised">
##META=<ID=SwitchedAlleles,Number=1,Type=Integer,Description="Total number of variants strand switched">
##META=<ID=TotalControls,Number=1,Type=Integer,Description="Total number of controls in the association study">
##META=<ID=TotalCases,Number=1,Type=Integer,Description="Total number of cases in the association study">
##META=<ID=StudyType,Number=1,Type=String,Description="Type of GWAS study [Continuous or CaseControl]">
##SAMPLE=<ID=ukb-d-20544_1,TotalVariants=9938650,VariantsNotRead=0,HarmonisedVariants=9938650,VariantsNotHarmonised=0,SwitchedAlleles=209,StudyType=Continuous>
##contig=<ID=1,length=249250621,assembly=HG19/GRCh37>
##contig=<ID=2,length=243199373,assembly=HG19/GRCh37>
...
##contig=<ID=22,length=51304566,assembly=HG19/GRCh37>
##contig=<ID=X,length=155270560,assembly=HG19/GRCh37>
##contig=<ID=Y,length=59373566,assembly=HG19/GRCh37>
##contig=<ID=MT,length=16569,assembly=HG19/GRCh37>
##contig=<ID=GL000207.1,length=4262,assembly=HG19/GRCh37>
##contig=<ID=GL000226.1,length=15008,assembly=HG19/GRCh37>
...
##contig=<ID=GL000192.1,length=547496,assembly=HG19/GRCh37>
##gwas_harmonisation_command=--json /mnt/storage/private/mrcieu/research/scratch/IGD/data/dev/ukb-d-import/processed/ukb-d-20544_1/ukb-d-20544_1_data.json --ref /mnt/storage/private/mrcieu/research/scratch/IGD/data/dev/QC/genomes/b37/human_g1k_v37.fasta; 1.1.1
##file_date=2019-11-25T15:12:05.258952
##bcftools_annotateVersion=1.9-74-g6af271c+htslib-1.9-64-g226b4a8
##bcftools_annotateCommand=annotate -a /mnt/storage/home/gh13047/mr-eve/vcf-reference-datasets/dbsnp/dbsnp.v153.b37.vcf.gz -c ID -o /mnt/storage/private/mrcieu/research/scratch/IGD/data/dev/ukb-d-import/processed/ukb-d-20544_1/ukb-d-20544_1.vcf.gz -O z /mnt/storage/private/mrcieu/research/scratch/IGD/data/dev/ukb-d-import/processed/ukb-d-20544_1/ukb-d-20544_1_data.vcf.gz; Date=Mon Nov 25 15:37:00 2019
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  ukb-d-20544_1
1       692794  rs530212009     CA      C       .       PASS    AF=0.111693     ES:SE:LP:AF:SS:ID       0.000379398:0.000798869:0.197332:0.111693:117706:1_692794_CA_C
1       693731  rs12238997      A       G       .       PASS    AF=0.116734     ES:SE:LP:AF:SS:ID       0.00058531:0.000756397:0.357493:0.116734:117706:rs12238997

For the first part of question 2:  if there are other recommended sources like dbGaP.
Please refer to this Databases of GWAS summary statistics:

Database
Content

GWAS Catalog
GWAS summary statistics and GWAS lead SNPs reported in GWAS papers

GeneAtlas
UK Biobank GWAS summary statistics

Pan UKBB
UK Biobank GWAS summary statistics

GWAS Atlas
Collection of publicly available GWAS summary statistics with follow-up in silico analysis

FinnGen results
GWAS summary statistics released from FinnGen, a project that collected biological samples from many sources in Finland

Pheweb.jp
GWAS summary statistics of Biobank Japan and cross-population meta-analyses

https://www.nature.com/articles/s43586-021-00056-9/tables/3

Thanks for pointing to that paper, very helpful +1

However, the column names won't match
exactly. I know that I can use BETA instead of OR, and FRQ might
correspond to MAF, but I can't seem to find the imputation information
score INFO.

Yes. There is not a standard format that everyone reports GWAS summary statistics in, unfortunately. First, check the paper that you are sourcing it from to see if they already filtered on imputation quality. Do ctrl+f and look in the methods section for imputation information (you may also try searching for imputation tools, like minimac or Beagle).
The good news is that if the sumstats have many alleles, imputation accuracy doesn't matter much. Errors get averaged out and each error only has a small effect since there are many alleles. Usually it is expected that reported sumstats will be pre-QC'd in some way, especially if they lack INFO. It would be safe to proceed without the INFO filter step.

The next step in the tutorial involves filtering based on specific
criteria, such as MAF and INFO values. However, I can't proceed
because these columns are missing.

Sumstat formatting is whatever the author wants, and can be a bit annoying to parse. However, FRQ is basically MAF.
See here (ctrl+f "frq"):
https://bioconductor.riken.jp/packages/3.14/bioc/vignettes/MungeSumstats/inst/doc/MungeSumstats.html
https://academic.oup.com/bioinformatics/article/37/23/4593/6380562

I would like to primarily use the terminal and only turn to R for
statistical analysis. Could you guide me on the correct approach to
convert the .vcf file to .txt so that I can continue following the
tutorial with PLINK?

Did you try unzipping and simply replacing .vcf with .txt in the filename?

General approach In addition to addressing the conversion issue, I'm also seeking guidance on the appropriate data sources for genetic
analysis. I'm curious about whether OpenGWAS is a reliable source for
publication or if there are other recommended sources like dbGaP.
Additionally, what is the ideal data format to obtain?

This aligns with your first question. If you use a database, sumstats will be standardized (usually) to the database's common format, so there will be less confusion over inconsistent column names. The GWAS catalog pretty good in my opinion, and all of them are acceptable.


text:
How to BLAST a protein against a large set of organisms?
I have a list of organisms with their full scientific names and Taxonomy IDs, total ~1500. I want to BLAST my protein against this 1500 genome or proteome. How can I do that?
Here a sample protein that I want to search the following database:
http://togodb.org/db/thermobase
Please note that database can be downloades as CSV file. I am very comfortable using web tools or programmatic tools such as Perl, Python, bash, or Entrez Programming Utilities.
>QCC44558.1 chitinase [Halobacterium salinarum] TaxID 641522631
MPHDRRSYLRTSSAVIASLLAASTPTSAADTPPEWDPDTVYTDGDKATFDGYVWEAKWWTKGDKPGADEW
GPWTQRRPVDGSPGGPTAAFTVTDTTVEPGTAVTVDASDTTGSVDSYEWAFGDGTTAAGVTASHTYDAAG
EYTIELTVTTGDGTTDTATQQITATTSPGDDEFKVVGYYPSWKGTDDYDFYPADVPFDQVTDVLYAFLDV
QPDGTVVLPDSDVDHESLLASFASLKQDRAADTRLKLSIGGWGLSPGFEDAAADQASRERFATTAVDLMR
TYDFDGIDVDWEHPGPRRGKCECGSAQGPANHVALLETVRDRLNDAEVEDGRTYDLSVANGGSDWNAALI
NHREVAAVVDDIYMMAYDFTGVWHGTAGLNAPIYGTPPDYPPSGDAQQYTLETTLTIWKEQGYWVDWMEW
EDHGDPVDDPGTLVLGMPFYGRGCNVENGIWDTFSLSEWQQGDPKYQNDVIPPGTWNDLRGPDDANTGAF
DYGDLAANYEGADAWTKQRNEQGGVPYLWNDSEGVFISYDDPTSIAAKVELAVAEDLGGVMIWELSQDYD
GTLLETINQHTP

@M__ Yes you are assessment is correct. I will consider breaking it down. I guess I have to build a local BLAST database first using TaxIDs. I will give try and maybe delete my question.

I will propose a different solution over the weekend. It's a lot easier than performing a local blast. My strong advice is to take the easy route. Whilst @terdon and I have done this style of calculation - there's much easier ways. This is because for major repeat projects - this level of coding is worth while - for what appears a one-off blast its not worth it. Also interogating NCBI via its link system is never fun (it's not standardised, i.e. thier databases don't use universal foreign keys). Basically solutions around remote blast only - there's a lot of caveats here and I'll describe them

My solution is simply.

perform a web-based Blastp using the following parameters:

most dissimilar
5000 hits (maximum limit)
Restrict to Archaea (which is what you samples appear to be)*

If you hit the limit of 5000 results, I would repeat the search for each class or phylum within Archaea. I can't remember if it was re-classified as a kingdom or was just another phylum. Concatenate all hits into a single fasta file.

Reg-ex your extremophiles database (personally I'd use pandas), so that the species name of each bacterium is recovered no other information is present. Or this is an easily accessible value within a hash/database (this is a better solution).

Use Biopython to loop through your blast fasta file and filter using re search against the species from your database. If you get a hit save that sequence.

Ideally, the point 2 and 3 would be a pandas merge operation, but looping is easier.
*, If there are extremophiles which are not Archaea (I don't know much about them), then take each general taxonomic group and repeat the operation for the search criteria "restrict the blastp to ... specific group".

Issues The issues with command line blastp, not via the web, is that access is greatly restricted during US working hours. If you are searching more than 100 hits and carrying out the search during the working day US time - they will not let it through. 100 hits or less is fine however. For >100 hits (e.g. 5000) it's got to be done well outside US working hours/ weekends. Web-based blastp doesn't carry this restriction of access during US working hours. Basically, the API route is deprioritised for the web-based approach. I assume this is to stop bots - albeit of course a bot can automate through NCBI's webpage.
Note The original solution by @terdon might easily work as a remote blastp because each query is restricted to <100 - so it will always get through during working hours. What I am saying is if you're using Blast API it needs a carefully defined strategy, else shift to the web for blanket approaches.
Entrez operations from a commandline do work ok during US working day time - albeit it's got to be coded to deal with "random" server refusals otherwise it will break the code (sporadic queries are refused). This is just my recent experience nothing more.
Just to be clear - this is potentially a medium sized project. Once the strategy has been decided then specific code that supports this approach can be given in separate questions.

I deleted my answer because I couldn't get the entrez query to work well. Too many of the species names were problematic.

@M__ your comment on the Issues is revelation to me as in the past I tried to BLAST  >100 proteins at the same time. And sometime it worked other time it did not. Now it makes perfect sense. I have also submitted a ticket to BLAST team, if they offer another solution I will share.

NCBI Entrez not working @terdon ... I've encountered bizarre non-universal foreign "link" key issues (I assume they're using SQL). In bacteria this is terrible: there were separate data and link structures present for one protein database and not the other, particularly how they linked to nucleotides ... it's disjointed. As a user I had to write two separate solutions in the code (I wouldn't mind but they never tell you either, got to figure it out). Whats happened is separate teams worked on "biologically the same" problem. The result is a committee meeting.

@Supertech ticket to blast support I tried this too. I hope you get better success than I had.  You said "[my] comments ... revelation", thanks, but upvotes and "accepts" are much more welcome than "revelations".

In the interface of blastp, you can limit the search space to organisms you want to. Yo can do this by introducing name species or Tax IDs, you can add more than one, but I bet you don't wanna introduce the 1,500 one by one.
The only other way to do this, is by using blastp in its command-line version. You would need to build a custom database to be used as the searching space, and then run blast on your computer or server using this database. Bulding it is a pain in the ass, as far as I know. You need first to download all the proteins in all organisms in fasta format and then use this compilation as a database (makeblastdb command). I am not aware of a way to download all proteins from various organism given a list of taxonomic IDs, but if any, it has to be NCBI, Uniprot or similar. Then you can use blastp locally together with your custom database.
But why limiting the search space is such a way? You are biasing your results to what you want to get, instead of finding the stuff that actually is closest to your protein of interest. My suggestion is to use the online version of blastp with one of the various databases available and rather than pre-selecting organisms, you can exclude the ones from which you definitely don't want to get results from, and you can also exclude entire families/genera such as primates (taxid:9443), to name an example.

There are many valid reasons why you would want to limit the results. For instance, if you are looking for the protein's homologs in particular taxa, you don't care about what is closest.

If you have noticed, database I want to use contains only (hyper)thermophilic organisms. That is exactly what I want.

There should be a way of downloading genomes based on the TaxIDs. I will send a ticket to the NCBI support team. In the past they have been very responsive. Maybe a feature they want to add if is not already possible.

@terdon, True, there are such cases or others. I still would suggest to first try to play with the "add/remove organism" boxes in the interface and then parse the blast output in order to search for the organisms of interest.


text:
Q-Q plots and the genomic inflation factor (k)
I have done the limma association analysis using below model design but not sure how to  evaluate Potential systematic biases using Q-Q plots and the genomic inflation factor (k). Could you please help with this?
#model matrix
var<-model.matrix(~BMI + as.factor(Gender) + Age +CD8T +CD4T +NK + Bcell +Mono ,data=targets2)
fit<-lmFit(mval,var)
fit2<-eBayes(fit,trend=TRUE, robust=TRUE)
probe<-topTable(fit2,adjust="BH",coef=2,num=Inf)

write.table(probe, file = "limma.output.txt")
data <- read.delim("limma.output.txt")
head(data)

               logFC       AveExpr.    t         P.Value  adj.P.Val         B
cg03546163  0.41018776  -0.8013218  7.85104192  1.13E-12    8.65E-07    12.3051291
cg17944885  -0.4884859  2.31974279  -5.7319645  6.33E-08    0.02431349  5.15237109

#For Genomic inflation factor

# For p-values, calculate chi-squared statistic
chisq <- qchisq(1-data$P.Value,1)

#Calculate lambda gc (λgc)
median(chisq)/qchisq(0.5,1)

#QQ PLOT

library(gaston)
pdf("qq.plot.pdf", width = 12, height = 12)
qqplot.pvalues(data$P.Value)
dev.off()

What have you already tried? see e.g. here and here that i found by google quickly.

thank you so much @MaximilianPress, I am not sure do we need to filter limma output based on adj.P.val before running QQ plot or not? and second I am not sure how to proceed after this step.

adj.P.val computes q-values (or applies some other correction), so its output are not really p-values. You will need to plot the p-values themselves, not a transformation.

thanks@MaximilianPress, I updated code above for QQ plot and genomic inflation factor. I hope its correct?

I don't see any obvious issues at this point but I'll let other possibly more expert people comment.

You do not need to filter the limma output based on adj.P.Val before generating the Q-Q plot / getting the inflation factor. Your code looks good to me. Is there anything unusual about the Q-Q plot or inflation factor that makes you think something went wrong?

Thanks for all help.


text:
Is it possible to, all in memory without writing a file, given a list of AlignedSegments and a bam header, build a bam, index and get coverage?
My current job is rate-limited chiefly by the IO overhead of writing reads to their corresponding bam files (they are split up by both cell barcode and contig). The only reason I need to write these mini-split bams is so that I can then use the pysam coverage function to assess read depth at certain locations. I already have all of these reads in memory as a list of AlignedSegments. Is it necessary to write them to a .bam file in order to be able to index and use the coverage function? I feel like there should be a way to have the AlignedSegments be treated as a bam object in memory, possibly even indexed in memory, and then treated as if they were just streamed from a bamfile with an index in order to be able to use the coverage method. Thanks in advance, if this is possible it would increase my processing speed x10. Or, given a list of AlignedSegments, is there a way to get coverage at specific positions across all those AlignedSegments in my list (that would get me what I need).

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.

What operating system are you using?

My operating system is linux

Assuming you are on an operating system that supports them, like Linux, you can try using named pipes in a ramdisk or tmpfs. Not the most elegant approach, but I don't know if what you are asking for is possible (it may very well be, I just don't know) and this can be a decent workaround.
First, check that you do have a mounted tmpfs. On my system, I have a few:
$ mount | awk '$1=="tmpfs"'
tmpfs on /dev/shm type tmpfs (rw,nosuid,nodev,inode64)
tmpfs on /tmp type tmpfs (rw,nosuid,nodev,size=16278892k,nr_inodes=1048576,inode64)
tmpfs on /run/user/1000 type tmpfs (rw,nosuid,nodev,relatime,size=3255776k,nr_inodes=813944,mode=700,uid=1000,gid=1000,inode64)

These filesystems are actually stored in RAM and should be significantly faster. For instance, on my machine, which has an SSD as the main storage, I still get significant speedup when using /dev/shm:
## SSD
$ time seq 1000000000 > file
real    0m13.687s
user    0m5.623s
sys     0m7.516s

$ time wc -l < file
1000000000

real    0m4.972s
user    0m0.101s
sys     0m2.248s

## tmpfs, RAM
$ time seq 1000000000 > /dev/shm/file

real    0m6.755s
user    0m4.982s
sys     0m1.763s

$ time wc -l < /dev/shm/file
1000000000

real    0m1.024s
user    0m0.120s
sys     0m0.902s

As you can see, using /dev/shm is significantly faster. And using named pipes can give you a little boost too:
$ time ( seq 1000000000 > /dev/shm/file;  wc -l < /dev/shm/file)
1000000000

real    0m8.341s
user    0m5.277s
sys     0m3.044s

$ rm /dev/shm/file && mkfifo /dev/shm/file
$ time ( seq 1000000000 > /dev/shm/file &  wc -l < /dev/shm/file)
1000000000

real    0m6.564s
user    0m5.426s
sys     0m3.965s

So, this is essentially a way of doing everything in RAM without needing to code it as doing everything in RAM. You can create these files (the FIFOs, the named pipes) and treat them as files, but since you're creating them in a tmpfs filesystem, they're actually in RAM and by using FIFOs you not only get a further speed jump but also minimise the actual data that need to be stored.

This seems like it gets at what I'm trying to accomplish. I will try it! Thank you.


text:
Making a DNAStringSet object
I wish to download a fasta sequence from NCBI using "rentrez" package and then save it in an object, then I want to have this object as a DNAStringSet object to work on the sequence using "Biostrings" capabilities.
I leave my incorrect code here so I would be grateful if you could help me out.
cov2 <- entrez_fetch(db="nucleotide" , rettype = "fasta", id = "NC_045512")
cov2 <- readDNAStringSet(cov2)
...
Error in .Call2("new_input_filexp", filepath, PACKAGE = "XVector") :
  cannot open file '>NC_045512.2

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.

Check what the value of cov2 is after the first command. entrez_fetch() returns the fasta record, but readDNAStringSet() requires a file path as input.

Your issue is that entrez_fetch() returns the fasta record(s) as a single string with the definition lines and sequences concatenated. Also, newline characters (\n) appear in place of line breaks, with double newlines (\n\n) separating different records. Here's a function that creates a DNAStringSet given a list of nucleotide identifiers without writing each fasta to a file:
entrez2dss <-
  function(id_list) {

    require(rentrez)    # v1.2.3
    require(Biostrings) # v2.66.0
    require(stringr)    # v1.5.0

    # fetch all sequences, trim empty elements
    raw_ <-
      entrez_fetch(db = "nucleotide",
                   rettype = "fasta",
                   id = id_list) |>
      str_split_1(pattern = "\n\n") |>
      head(-1)

    # get sequence names
    names_ <-
      gsub(pattern = "\n.*",
           replacement = "",
           x = raw_) |>
      gsub(pattern = "^>",
           replacement  = "")

    # get sequences
    seqs_ <-
      sub(pattern = ".*?\n",
          replacement = "",
          x = raw_) |>
      gsub(pattern = "\n",
           replacement = "")

    # create named DNAStringSet object
    dss <- DNAStringSet(x = seqs_)
    names(dss) <- names_
    dss

  }

seqlist <- c("NC_045512", "OR558592")
cov2 <- entrez2dss(seqlist)

The previous error message was because the entrez_fetch() function returns the sequence as a character string, not a file path. To save the sequence to a file, I can use the writeLines() function. For example, the following code will save the sequence to a file named cov2.fasta:
# Download the FASTA sequence
cov2 <- entrez_fetch(db = "nucleotide", rettype = "fasta", id = c("NC_045512", "OR558592"))
# Save the sequence to a file
writeLines(cov2, "cov2.fasta")
# Read the FASTA sequence into a DNAStringSet object
cov2_dna <- readDNAStringSet("cov2.fasta")


text:
Asterisk (*) calculation method
Does anyone here know how to calculate a value for the asterisk (*) code that appears in substitution matrices?
From my observation, to all pairs with a one asterisk, the lowest value from the matrix is set. For ** it's 1. Example: BLOSUM62.
But, it's just my observation, I couldn't find any source or explanation how it was calculated (for the exemplary BLOSUM62 or any other substitution matrix).
I'd appreciate sources, especially scientific papers.

Can you please explain what you mean by "a value"? You already gave the values: 1 for another stop, and -4 for any aminoacid in the example you linked to. What else do you need?

I want to know what's the formula / method behind these values / scores.

OK, in which matrix exactly, then? PAM? BLOSUM? Which PAM? Which BLOSUM? do you know that the asterisk stands for a stop codon?

Ye, I heard about it.

BLOSUM would be the best, but as I couldn't find it for a long time, so I would appreciate any source.

Its really easy. No formula is needed.
Protein matrices are based on observed frequencies from nucleotide alignments for a variety of genes. They are not based on a priori formulas (it's not like nucleotides). Later matrices incorporated phylogenetic criteria using maximum likelihood - it's really complicated how it works.
Blosum used odds-ratio (OR) I think without a phylogenetic tree. I suspect raw OR. Explaining with/without a phylogenetic tree is a full lecture in phylogenetic theory, just think it'll be a raw value calculated directly.
So lets just run with raw OR based matrix ...
So ... * is a stop-codon it carries the lowest weight in the matrix because each protein only has one of them. So a stop-codon frequency in comparison to any amino acid is going to be extremely low - because for that aligned position there's only stop-codons (i.e. no other amino acids).  Thus its -4 here (lowest value in the matrix).
What about the 1?
However, when you compare a stop codon frequency ... between proteins, well its occurrence is always 1 to 1, right? Each protein has one stop codon, compare it to another protein ... thats got one stop-codon too. So one million proteins which are nicely aligned to be homologues (thats how its calculated) have one million stop codons ... thus its 1 no matter how many protein there are in the alignment. The alignment position is always the same because stop codons (at least in Blosum) are always homologous.
Once more, comparing a stop codon (in a protein), with the other stop codon in another protein some the same homologous "amino acid" site - no matter which protein it is - its always going to be 1 because its a universal feature of all proteins.
Thats how it works and I would simply cite

Henikoff, S.; Henikoff, J.G. (1992). "Amino Acid Substitution Matrices from Protein Blocks". PNAS. 89 (22): 10915–10919

I've checked the paper and the authors didn't use stop codons, nor did they use any phylogenetic criteria (that was easy to guess). I know they didn't include stop codons because because,

The Sigma function for matrix operations has a 20 above it (i.e. 20 amino acids), there's no allowance for stop codons.
The matrices they present omit stop codons.

When NCBI are leveraging this matrix they need answers for all possibilities, the idea that X (any amino acid) or * is omitted would not be cool for them. Thus they could have approached the authors to fill in this information, this would probably mean redoing the alignments to include the stop codon (which the authors wouldn't have liked), or NCBI could have added it after the event (its a pretty safe addition) because the log-ratio of an invariant site is always 1. The authors/NCBI would just need to be sure the  stop codon was in the same position for all protein alignments in Blosum.
You might find the wiki more understandable  https://en.wikipedia.org/wiki/BLOSUM

Thank you for this detailed answer. I really appreciate it!

Would you be so kind to point out a paragraph containing this information? I know this paper and couldn't find it. Or just a definition that says an asterisk is a stop-codon, because the rest can be inferred out of it, I assume.

Thanks @maciejwww I answered you question below the '-------------'

Thank you so much!


text:
How to use HGTECTOR tool for identifying HGT genes in a bacterial genome?
I wanted to know the HGT genes in a Metabacillus strain and I noticed that several papers use HGTector for this purpose. I installed HGTector through conda and tried to download the database for Bacillota using the example codes given on the hgtector github page -
hgtector database -c bacteria -o db1 -t 1239

The code for Bacillota is 1239. It took an awfully long time (approx. 13h) for the download to end and several warning messages were displayed saying so and so files could not be downloaded. When I checked "db1" only 1000 or so files were there, signifying that the download was incomplete as there are 9112 files under Bacillota. I am not sure if the problem is with my code or because there is not enough disk space or something else. Following are my system and Wifi details :
Lenovo Ideapad, 16GB memory
12th Gen i5 - 1235U x 12 Processor
1 TB Disk capacity
OS = Ubuntu 22.04.3 LTS
Wifi speed = 32Mb/s

I have posted this on github about 4 days ago but no response yet. If  my disk space is not enough could anyone suggest any other tool to study horizontally transferred genes? I have scanned the web but did not find anything satisfactory.
I ran df -h and it shows ~700 GB of free disk space even after the download of 1000 files which should be around 2 GB.
I've tried again at 'family' level with similar results. I don't think I can use HGTector for now. What do you think of RecentHGT? This is the paper for it - link
Thank you

Have you tried running the download in a different way? Like, can you download the db directly and not through the tool? They provide links on the github for prebuilt databases.

Yes I did, but I think the disk space is not enough. My system does not show any error related to this but either way, the download stops midway with not so meaningful  warning messages. So I just gave up.


text:
how to add a custom elements into PyMOL visualization?
I would like to add some 3D elements (balls and lines) to the PyMOL visualization to highlight some "holes" in the structure. Say I fetch 1tqn and I want to add a yellow ball somewhere (to specific coordinates, say [-19, -23, -14]) in the view. Or a bar from - to coordinates.

I was going through wiki, did a online search but have not found a way how to add these simple 3D entities into the scene. I guess it should be trivial and I'm just using wrong queries. Because I need to add lots of these, I need to be able to do that via scripting interface.
Thank you in advance for any help!

see here https://stackoverflow.com/questions/2060582/draw-a-colored-sphere-from-cartesian-coordinates-in-pymol/76749167#76749167 Draw a colored sphere from cartesian coordinates in pymol  and here https://pymolwiki.org/index.php/Load_CGO and https://pymolwiki.org/index.php/CGO_Shapes

my attempt :
import pymol

from pymol.cgo import *

from pymol import (

                    cmd ,
                    )

print('########## PYMOL VERSION ##########################################')
print('         ',  cmd.get_version() )
print('###################################################################')

pymol.finish_launching()

cmd.fetch('1tqn')

cmd.remove('solvent')

pos = cmd.get_position()

r = 3

print('###############################################################################\n')

print('### POSITION : ' , pos,'  ###\n')

print('###############################################################################\n')

spherelist = [
    [COLOR,    1.000,    1.000,    0.000,
    SPHERE,   pos[0] , pos[1] , pos[2]  , r ],

    ]

for i in spherelist:

    a = str(spherelist.index(i))

    cmd.load_cgo(i, 'sphere_'+a )

cmd.zoom()

cmd.set('cartoon_transparency' , '0.5')

result:

Aaah, load_cgo() https://pymolwiki.org/index.php/Load_CGO ("Compiled Graphics Objects"). There are also other shapes: https://pymolwiki.org/index.php/CGO_Shapes. Thank you!


text:
How can I export a pruned phylogenetic tree in nexus format in R?
I have a phylogenetic tree with a lot of different populations and I wanted to removes most of them and apply new names to the labels. I wrote this code:
library(ape)
library(phylotools)
tree<- read.tree("Consensus_COMBINED_Cod3genesRatesC_trees.trees")

tree<- drop.tip(tree, tip = c("CeuMist2","CeuMist1","Ceu_PAN_2","Ceu_PAN_1","Ceu_PAN_3","Ceu_PAN_4"))
tree<- drop.tip(tree, tip= c("CcNav1","CcNav2","Cc_Mal_1","Cc_KalR2_2","CcRod1","CcRod2","CcRod3","CcRod4","CcHr1"))
tree<- drop.tip(tree, tip= c("CeST6_4","CeST6_5","CeAA7_1","CeST6_3","CeST6_2"))
tree<- drop.tip(tree, tip= c("Cg_MelR3_1","Cg_MelR6_2","Cg_MelR7_2","Cg_MelR6_1","Cg_MelR7_3","Cg_MelR5_1","Cg_MelR4_2","Cg_MelR1_1","Cg_MelR2_1","Cg_MelR4_1","Cg_MelR7_4","Cg_MelR5_2"))
tree<- drop.tip(tree, tip= c("ChTM5_1","ChTM4_3","ChTM4_1","Ch_MEN_5","Ch_MEN_13","Ch_MEN_8","Ch_LEV_3","Ch_LEV_5","Ch_LEV_4","Ch_MEN_6","Ch_MEN_10","Ch_LEV_2","Ch_MEN_4","Ch_MEN_14","Ch_LEV_1","Ch_MEN_11","Ch_MEN_3","Ch_MEN_9"))
tree<- drop.tip(tree, tip= c("Ci_MSP_2","Ci_MSP_1","Ci_DIA_1","Ci_DIA_2","CiPLAN1","CiPLAN3","CiPLAN2","Ci_MEN_1","Ci_MEN_2","Ci_MEN_12"))
tree<- drop.tip(tree, tip= c("CpGravR3","CpGravR4","CpGio_2","CpGravR5_2","CpGravR2","CpGravR7","CpGio_1","CpGio_3","CpGravR1","CpGravR6_2","CpGravR5_1","CpThe_1","CpSka_1"))
tree<- drop.tip(tree, tip= c("Cn_IG_2","Cn_IG_3","Cn_IG_1","LevantinaFC4779","Assyriella_naegelei"))

newlabels<- c("C. eucineta", "C. codringtonii", "C. elisabethae","C. gittenbergeri", "C. helenae", "C. intusplicata", "C. parnassia")
tree$tip.label<- newlabels

write.nexus(tree, "CodTree.trees")

When I use the write.nexus command I get the following error:
Error in y$tip.label : $ operator is invalid for atomic vectors

How can I circumvent this error and export the tree in nexus format?
I will also provide the tree here: https://drive.google.com/file/d/1lYmzW9htwaclC_yYs3fQom-Z3zbjk8wH/view?usp=drive_link

Hi @Nickmofoe I'll look at this in detail, but will not be until the end of the week (really busy week).

So after digging I found the solution. What worked for me was adding the file =  in the write.nexus command. The final line should be
write.nexus(tree, file =  "Cod.trees")

Ahhh! Right of course! Many thanks for posting back. Could you very kindly mark you post as "accepted". Admittedly, you don't any bonus rep, but it helps site statistics and it is definitely the solution. I completely missed the file =.


text:
GWAS in python?
A question similar to this has also been asked on Biostars
I understand GWAS is done with R, but are any written in Python?
I'm used to python and java, acclimating to R will take a bit of time and my advisor wants some data analyzed. I plan on using R but given the deadline, I want to finish the analysis in Python.

Please [edit] your question and give us (a lot, a huge lot) more context. Tell us what input data you have, what kind of analysis you want to do (do you have the results of the GWAS study or are you trying to run the study itself?), what species and anything else that could be relevant. Of course you can do it in python, you can do anything in any programming language, so you really need to narrow this down to something very specific that can be answered in a few lines.

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.

I would adapt PLINK. Whilst Biopython pop gen could be used it is better to use a standalone package.

You can use Python to wrap around command line tools

It doesn't (really) matter what language a piece of GWAS software is written in. SAIGE is written in a mix of R and C++. BOLT is written in C++. You shouldn't be writing a GWAS analysis from scratch, as it ends up being a very complex piece of software. You will almost certainly be interacting with the command line, when running a GWAS, so my advice is to find the software which is most suitable for your piece of analysis, not one which is written in a programming language you like.


text:
1.What size PCR product will be generated using the primers (bold and underlined) in the sequence below?
I would like to know the steps of this question and how to determine the forward and the reverse primers and also the length of the sequence, without using NCBI. If you could explain and give an answer step by step I would appreciate it. Thank you

Wouldnt it be a question for chemistry or biology stack ?

You are the second person to have raised this @pippo1980 exact point, Biology SE was suggested. If you flag this formally it helps to decide about migration.

I think it should be on SE Biology. Trying to flag with this.

This question tells me that you do not understand how PCR works. You should figure it out. Use Google search.

Just to note, the question has now received 2 closed requests.

Without giving too much of the "what they want you to do yourself" bits of this away, here are a few explanations of the image that may help you answer this question:

As mentioned in the text, the forward and reverse primers are underlined in the image.
The - and + symbols below the forward-direction sequence help determine the positions of the bases above and below. A + appears every 10 bases.
The first base in the sequence is at position 1441 (as indicated by the numbers) to the left of the sequence. The first + to the right of that first base represents position 1450.
Primers are usually represented as a sequence in a "5-prime to 3-prime direction", which is almost always the left to right direction of the top sequence in an image like this, and the right to left direction of the bottom sequence.
The length of a PCR product includes the primer sequence.

For help beyond that explanation... this looks like a canned question. Unless you give more context around what bioinformatics project you're working on and how this question relates to your project (i.e. explain why it's a real-world problem rather than a text-book question), I'm going to leave answers to the specific questions asked as an exercise to the reader.

So in this case what would be the length?

I'm happy to help you understand the question, but I don't want to help you answer this question beyond that understanding. Other people are welcome to do that if they desire.

Thank you for that! But as the questions says, I need to determine the size of the PCR product. Is there a formula for that?

@Teo you should be able to figure out the formula from what gringer has told you. We are not a homework writing service, our job is to help you understand so you can figure it out. So please read through what gringer has written, and see if you can figure it out. Then, if you still need help, update your question and add what you think the solution is and how you found it and then we can help you take it from there. If the information given here isn't enough for you, then update your question and add what you understand and what you feel is missing.

something about 1.5MDa in weight ?

@Teo Indeed there is a formula: b - a, where a is the start and b is the end point :)

@Teo if I start reading a book on page 20, and then close the book after page 30 and go have a cup of coffee, how many pages have I read?


text:
How can I incorporate wildcards in Snakemake in an R script?
I am facing the following issue:
I have a rule in Snakemake that looks something like this:
rule somerule:
input:
    tables = expand("results/{tables}_table.txt", tables = ["1", "2"]),
output:
    edit= expand("results/{tables}_edited.txt", tables = ["1", "2"]),
conda:
    "../envs/r.yaml"
script:
    "../scripts/script.R"

The script used in this rule is as follows:
library()

tab <- read.table(snakemake@input[["tables"]], sep = "\t", header = TRUE)
edit <- some_function(tab)
write.table(edit, file = snakemake@output[["edit"]], sep = "\t", quote = FALSE, row.names = TRUE, col.names = TRUE)

When I execute the code, I get this error:
Error in file(file, 'rt'): invalid 'description' argument
Called from: read.table -> file
Execution halted

However, when I input the data individually (not as wildcards), the script works fine and produces the correct output files. Where could the problem be? My suspicion is that Snakemake is passing all the input data to the R script simultaneously instead of running them individually, but I'm not sure if that's the case and even less sure how to fix it.

Your problem is that Snakemake is passing all the input data to the R script simultaneously instead of running them individually :)
As written your rule takes a list of input files and gives a list of output files (expand is just a formatting helper function that returns a list based on its arguments).  Instead maybe have one rule that calls the R script with one input and one output, and another rule that asks for all those outputs as an easy way run the whole thing at once.
Maybe something like this?  I changed tables to table here just to emphasize it's one by one.
rule all:
    input: expand("results/{table}_edited.txt", table = ["1", "2"])

rule somerule:
    input:
        table = "results/{table}_table.txt"
    output:
        edit = "results/{table}_edited.txt"
    conda: "../envs/r.yaml"
    script: "../scripts/script.R"

thank you so much. That worked!


text:
Does pooled sample genotyping microarray analysis exist?
I want to genotype samples in small batches with a common microarray (e.g., the GSA). However, the cost is very high for small batches.
With a higher sample throughput/volume, per sample cost becomes more reasonable as the processing cost for a single run is spread out over more samples.
I'd like to process no more than 16 samples per run, ideally even one or four samples per run. (Being able to submit a single sample for microarray analysis without having to pay to run the rest of the empty chip would be great.)
It is possible that some company or core offers "pooled" microarray services, where I could submit a small number of samples that would be processed with other groups' samples on the same microarray run. This would mean there are more samples per run, so it could potentially cost much less.
Looking on core and company websites, I don't see this offered anywhere.
Does this service exist?

I understand the question entirely. My personal view is this is primarily a wet-lab question. Companies will do mass chip printing to reduce costs, which are stored. Grouping between independent clients I would suggest is avoided because of the risk of any form of "spill-over". For example, chips get mislabelled and you get someone else's results.  For the company it's a disaster and they lose client confidence. Better to work alongside colleagues with similar aims as you for the "pool".

The experimental problem with microarrays/chip technology is the variation in quality. Statistically this can be corrected, but wet-lab perspective is they don't like any inclination of sub-optimal quality.

@M__ Yes, that makes sense. I imagine that's why most cores don't operate like that

This can also be done at the sender (rather than the chip runner) level. Multiple senders get together to send samples


text:
Cannot obtain alignment summary after running Bowtie2
I am aligning my Small RNA Seq data with Bowtie 2. Although the alignment performs well, the only information I obtain after finishing running the alignment is the following:
## Time loading reference: 00:00:00
## Time loading forward index: 00:00:00
## Time loading mirror index: 00:00:00
## Multiseed full-index search: 00:00:00
## Time searching: 00:00:00
## Overall time: 00:00:00
## Time loading reference: 00:00:00
## Time loading forward index: 00:00:00
## Time loading mirror index: 00:00:00
## Multiseed full-index search: 00:00:01
## Time searching: 00:00:01
## Overall time: 00:00:01
## Time loading reference: 00:00:00
## Time loading forward index: 00:00:00
## Time loading mirror index: 00:00:00
## Multiseed full-index search: 00:00:00
## Time searching: 00:00:00
## Overall time: 00:00:00
## Time loading reference: 00:00:00
## Time loading forward index: 00:00:00
## Time loading mirror index: 00:00:00
## Multiseed full-index search: 00:00:00
## Time searching: 00:00:00
## Overall time: 00:00:00
## Time loading reference: 00:00:00
## Time loading forward index: 00:00:00
## Time loading mirror index: 00:00:00
## Multiseed full-index search: 00:00:00
## Time searching: 00:00:00
## Overall time: 00:00:00
## Time loading reference: 00:00:00
## Time loading forward index: 00:00:00
## Time loading mirror index: 00:00:00
## Multiseed full-index search: 00:00:00
## Time searching: 00:00:00
## Overall time: 00:00:00

This is my command line to build the index:
bowtie2-build -o 3 -t 5 Nuevofastaorina human_mirs
And this is my command line to do the alignment:
bowtie2 -x human_mirs -U 1.lane.clean.gz -N 1 --norc --no-unal --no-head --ignore-quals –-quiet > sr1prueba.sam

Does anyone know how to obtain an alignment summary instead of the above information? What I mean by alignment summary is the following:
20000 reads; of these:
  20000 (100.00%) were unpaired; of these:
    1247 (6.24%) aligned 0 times
    18739 (93.69%) aligned exactly 1 time
    14 (0.07%) aligned >1 times
93.77% overall alignment rate

Please show the first few lines of your fastq input data (i.e. lane.clean.gz). This should help to diagnose the problem.

Remove the --quiet option, and capture STDERR as well as STDOUT:
bowtie2 -x human_mirs -U 1.lane.clean.gz -N 1 --norc --no-unal --no-head --ignore-quals 1>sr1prueba.sam 2>sr1prueba.summary.log

See also:
bowtie2 manual: alignment summary

Get samtools.
Then run:
samtools stats file.bam

and/or
samtools flagstat file.bam

For example outputs see here.


text:
Convert genotyping results to vcf or other plink formats
I used a Genotyping software called T1K to generate KIR genotyping, the output file contains genes and alleles for each sample, how can i convert the results to a vcf file or other formats which can be used in plink?
The file i have for now looks like this:
Gene_name Allele_1 Allele_2
KIR2DL1 KIR2DL1*004 KIR2DL1*042
KIR2DL4 KIR2DL4*005 KIR2DL4*005
KIR3DL1 KIR3DL1*005 KIR3DL1*005
KIR3DP1 KIR3DP1*006 KIR3DP1*010

My thought is to first convert the above file to a vcf file and use plink to generate bedbimfam from it, a standard vcf is like this:
##fileformat=VCFv4.3
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/1000GenomesPilot-NCBI36.fasta
##contig=<ID=20,length=62435964,assembly=B36,md5=f126cdf8a6e0c7f379d618ff66beb2da,species="Homo sapiens",taxonomy=x>
##phasing=partial
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele">
##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP membership, build 129">
##INFO=<ID=H2,Number=0,Type=Flag,Description="HapMap2 membership">
##FILTER=<ID=q10,Description="Quality below 10">
##FILTER=<ID=s50,Description="Less than 50% of samples have data">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality">
#CHROM POS      ID         REF   ALT    QUAL  FILTER   INFO
FORMAT       NA00001         NA00002          NA00003
20     14370    rs6054257  G     A      29    PASS    NS=3;DP=14;AF=0.5;DB;H2           GT:GQ:DP:HQ  0|0:48:1:51,51  1|0:48:8:51,51   1/1:43:5:.,.
20     17330    .          T     A      3     q10     NS=3;DP=11;AF=0.017               GT:GQ:DP:HQ  0|0:49:3:58,50  0|1:3:5:65,3     0/0:41:3
20     1110696  rs6040355  A     G,T    67    PASS    NS=2;DP=10;AF=0.333,0.667;AA=T;DB GT:GQ:DP:HQ  1|2:21:6:23,27  2|1:2:0:18,2     2/2:35:4
20     1230237  .          T     .      47    PASS    NS=3;DP=13;AA=T                   GT:GQ:DP:HQ  0|0:54:7:56,60  0|0:48:4:51,51   0/0:61:2
20     1234567  microsat1  GTC   G,GTCT 50    PASS    NS=3;DP=9;AA=G                    GT:GQ:DP     0/1:35:4        0/2:17:2         1/1:40:3

Any tools can be used to realize it? For now i only have Gene names and allele names for each sample, How to obtain other genetic informtation a vcf file needed?

Hi! I am not so familiar with this output format. Could you explain what we are looking at in the output sample?

@winni2k Thank you for your reply. I have added more details to my question as you suggested. Could you please take a look at it again?

I am not at all familiar with KIR genotyping, but is there a way of converting things like KIR2DL1*042 into chromosome, position, ref and alt? Is there some database that lists what these KIR notations actually mean?

I would also recommend searching dbsnp for the kir mutations.

your output file does not seem to match what is described in T1K docs. Are there other files output? T1K is supposed to output a VCF even.

@MaximilianPress  Thank you for reply. The output i posted is based on the genotype ouptput, The vcf output only contains noval SNP, i dont know how to obtain those info for other varients. Any methods to solve this? Thanks~

Genotype output according to the docs looks like this: gene_name num_diff_alleles allele_1 abundance_1 quality_1 allele_2 abundance_2 quality_2. That doesn't match what you show. But regardless, if they really don't want to give you non-novel alleles there is likely a reason, I'd suggest contacting the authors on how you might do that. I glanced at the code but was not sure how to get non-novel, though this line might be relevant.

@MaximilianPress Thank you, i really appreciate ur advice.


text:
How to incorporate negative controls in DESeq2
We are doing a comparison between two outcomes (positive and negative). We could not have any positive controls as we do not have any "control" data to set as baseline, either from literature or otherwise. No spike-ins have been used. We have some samples which had no resolved outcome, hence we are skeptic of using those as Controls.
We have few negative control samples which contain just the blank culture media with no RNA. We want to remove these signals from the RNA-seq data, but I would like to know the most ideal way to perform this. I would assume just normalizing and subtracting counts wouldn't be the most ideal.
So here's what the coldata looks like -
Sample_ID   Outcome
S1  Positive
S2  Negative
S3  Unknown
S4  Negative
…   …

My previous post has been about the same experiment
Here are the counts for the top 5 genes with high counts in negative samples (RIF1, GOLGB1, C10orf67, MAP3K15, F8):
Geneid  S1  S2  S3  S4  S5  S6  S7  S8  S9  S10 S11 S12 S13 S14 S15 S16 S17 S18 S19 S20 S21
ENSG00000189283.11  2284    2986    4577    2447    1290    1882    1953    1366    2443    3932    2640    1078    491 2919    1248    1627    1805    2094    2229    2789    190
ENSG00000173230.17  1000    1219    1324    1485    722 1142    870 943 1390    1617    1093    564 86  924 778 1059    1144    1248    1169    1828    50
ENSG00000180815.15  641 723 1118    1758    529 865 683 878 1554    2421    1571    414 136 1076    599 924 1098    1199    1405    1327    54
ENSG00000179133.15  369 433 592 488 305 366 335 252 455 2119    1357    118 159 314 200 284 376 333 517 550 17
ENSG00000256637.9   429 491 523 550 250 410 406 359 552 438 339 176 38  357 296 433 517 514 631 779 24

We plan to subtract the median count for these genes from the count matrix of all samples. If the value dips to negative, they will be replaced by 0. Is there a better way to handle background from negative controls? After all these signals are present in all samples, so I'm inclined to believe these signals are noise from the media itself, rather than anything biologically significant from the cells.


text:
biomarRT conversion of symbols
I wanted to use the code to convert mouse to human symbol, but I have problem with the mirror, I've tried to establish another one, but continue to have the same problem.
musGenes <- c("Hmmr", "Tlx3", "Cpeb4")

# Basic function to convert mouse to human gene names
convertMouseGeneList <- function(x){
    require("biomaRt")
    human = useMart("ensembl", dataset = "hsapiens_gene_ensembl")
    mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl")
    genesV2 = getLDS(attributes = c("mgi_symbol"), filters = "mgi_symbol", values = x , mart = mouse, attributesL = c("hgnc_symbol"), martL = human, uniqueRows=T)
    humanx <- unique(genesV2[, 2])

    # Print the first 6 genes found to the screen
    print(head(humanx))
    return(humanx)
}

Show us the error message.

My background in coding is more extensive in Python.  If you have Python skills, I can say for sure, use the Bioservices library:
from bioservices import BioMart

Install Bioservices via pip:
pip install bioservices

You can also use the web-based tool on Ensembl called BioMart.  It does the same thing, just easier to use if you continue to have difficulties with Python or in your case, R.  I looked up the R package for BioMart, it is called: biomaRt
Here's the link for the Ensembl BioMart web-based tool if you need to use it and are in a crunch for time on your project:
https://useast.ensembl.org/info/data/biomart/how_to_use_biomart.html#introduction
I hope this helps! :)


text:
Mitochondrial genes - TPM calculation bulk RNA-Seq
This question was also asked on Biostars
I was wondering if any of you have encountered a situation for bulk RNA-Seq where, possibly due to low sample quality or presence of dead cells, mitochondrial genes are expressed to a very large degree relative to other genes, thus skewing TPM values of all nuclear genes (by effectively scaling them down).
In such a situation, what could I potentially do to alleviate this beyond preparing fresh samples? Could I, for example, exclude 'MT-' genes from samples and then recalculate TPM based on this filtered set of genes?
My aim was to get some idea of within-samples expression levels, hence why going for TPMs as I know that counts normalised by DESeq2, edgeR and the like do not correct for gene length and are inappropriate for within-samples comparisons.
Many thanks in advance for any insight

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.

Yes, this is a common situation. Not just for mitochondrial genes, but for any gene that may be an outlier (e.g. the B-cell receptor gene in a contaminant/unexpected B-cell population). The way to get around this is to not use TPM for differential expression analysis.
DESeq2 by default uses a quantile-based normalisation method for bulk differential expression calculations, which works well when the majority of genes are not differentially expressed. When this assumption is broken (e.g. for a targeted panel), a control gene list can be supplied, or custom size factors can be applied to each sample.
More information about using DESeq2 for differential expression can be found in the DESeq2 manual:
https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
For comparing expression from different genes, I normalise variance-stabilised data using the longest annotated gene length from transcripts.
See this answer and this paper for more details.

Hi, thanks very much for your response. I should have explained this at the outset, but my aim was to get some idea of within-samples expression levels, hence why going for TPMs as I know that counts normalised by DESeq2, edgeR and the like do not correct for gene length and are inappropriate for within-samples comparisons.

Added a little more detail about how I do it, using the VST from DESeq2.

Thank you, appreciate the help I will take a look at those links


text:
Producing just network files from BiG-SCAPE
I was wondering if anyone knows off-hand if there is a way to run BiG-SCAPE without producing the svg or html files? All I want is the networking files (which I can play with in CytoScape).
There isn't anything obvious from reading the parameters section in the wiki but perhaps I just missed it? Maybe like a minimal mode or something similar?
Normally I would simply delete them, but I'm running around 20 000 BGCs through it. And I thought it would save time and disk space if those files are just never written out. I'm simultaneously trying out bigslice but the documentation is surprisingly tricky to follow.


text:
Sequence Alignment for sequences with the same length
I am doing research on a new method of optimizing sequence alignment process (Needleman - Wunsch algorithm) but the idea would only work with sequences that have the same length. I am wondering if that will have any real-life utility. Specifically, how often do we have those situations where we are aligning same-length sequences (or at least sub-sequences) and what are those situations?

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.

This is an open-ended question that will lead to opinion-based answers, which doesn't fit well with the Stack Exchange format. Could you please change your question to relate to a specific situation or problem? If that can't be done, I recommend instead asking this question on a more discussion-based site like reddit.

Why do you choose the Needleman-Wunsch algorithm?  Are you thinking of using a computer programing/ AI/ML approach to this question and project?  When you mention algorithm, I think more of a Python AI/ML approach.

I see ... the application is if two genes had the same number of indels* at different homologous positions within the gene. This happens but equally one gene might have an indel and another gene might lack it, in which case these genes would be different  lengths.
So just to reiterate your method would find genes with an identical number of indels at different homologous positions in a gene. Its application is protein modelling. This is because indels cause a large change in the entropy of a protein, so even if these two genes had identical sequence (barring the indels obviously) they would be structurally different.
Beyond this however the application is limited because apart from the known changes to protein structure, little is understood about quantifying the function of indels. Certainly in phylogenetics we specifically avoid them, because we struggle to fit them into a known model.

*, indel = insertion deletion event. In a protein they occur in triplets.

My 2 cents: Sequence alignment is a well studied problem in bioinformatics with some great existing methods that do handle uneven length sequences (i.e. indels). With tons of genomic data the trend these days is towards faster methods with more heuristics like pseudo-alignment or k-mer based approaches rather than theoretically optimal but slow and unscalable approaches like Needleman Wunsch. There might be some niche applications on DNA motifs or short peptide analyses but I think these are quite limited.


text:
getting NCBI TAXname by R
I have a data frame on R, called qIns8, with 14 columns and 204 rows. the column called "subject tax id" has NCBI taxIDs for the corresponding row, I need a script to add another column to the data frame called Organism for which the script has to get the organism name (NCBI taxname) from the column subject tax id.
I already tried this script using Taxize package, but it did not work:
getTAXname <- function(tax_id) {
  name <- tryCatch({
    tax_name <- sci_name(tax_id)
    return(tax_name)
  }, error = function(e) {
    return(NA)
  })
  return(name)
}

# Apply the function to each row of the dataframe
qIns8$Organism <- sapply(qIns8$subject_tax_id, getTAXname)

# Check the updated dataframe
head(qIns8)

I also tried this from Taxonomizr package and it did not work either:
getTaxonomy(qIns8$`subject tax id`, sqlFile = qIns8$`Organism`)

Please make sure to format code as code. See https://bioinformatics.stackexchange.com/help/formatting.

Please [edit] your question and explain what "didn't work" means. How did it fail? Was there an error? What error? Was there no error but you didn't get the expected result? What did you get? What were you expecting? Also, try to give us an example we can use. Simplify your script to the minimum code that reproduces the issue, provide us with a simple data frame that we can also load and try our solutions out with.

I developed the function and it worked perfectly,
Converst_taxID_to_taxName <- function(tax_id) {
  tax_info <- entrez_summary(db="taxonomy", id=tax_id)
  return(tax_info$scientificname)
}

# usage example
myDataFrame <- myDataFrame %>%
mutate(Organism = sapply(`subject tax id`, Converst_taxID_to_taxName))

# in case one tax ID has a problem or cannot be proceeded, it puts "NA" and automatically goes the next ID
Converst_taxID_to_taxName2 <- function(tax_id) {
  tryCatch({
    tax_info <- entrez_summary(db="taxonomy", id=tax_id)
    return(tax_info$scientificname)
  }, error = function(e) {
    return(NA)
  })
}


text:
Statistical approach to link DNA methylation with toxic element exposure and health outcome
I would be thankful to you if you can help with the statistical approach for case-control study to link DNA methylation (epic array) with toxic element exposure (arsenic) and health outcome (Cardiovascular diseases).
I do have data from 125 Cardiovascular diseases and 125 non-Cardiovascular diseases individuals matched with age, sex, and smoking.I do have phenotypes files in a csv format having sex, age, smoking, disease, samples id, iAs (arsenic) values in log format and then a separate files for methylation data.

Notes

The sample size is small.
Maybe structural equation modelling could be used to find CpG sites associated with arsenic exposure and disease condition as I have very small size samples
R is a preferred language.

This forum is not dedicated to general help requests and your post may be closed. Ask a precise bioinformatics question if you are encountering difficulties in your analysis, there is a lot of tutorials available on the internet.

What computer language are you going to utilize in this project?  Are you planning on an ML/AI approach?  That's a major question to ask yourself.  If you are planning on a statistical approach, I highly recommend R and then Python.  If your skills are advanced, Id try an AI/ML approach!  What databases/sources did you use to get the phenotype files in .csv format?  What data types are the csv values in?  Are they integers, decimals/floats, strings?  Should you consider a .tsv files format instead?  What's the purpose to your research question to arsenic and DNA methylation?

The question is part of a previous question here https://bioinformatics.stackexchange.com/questions/21823/mendelian-randomization . If they are taken together, plus comments it's answerable. I will make the appropriate edits shortly.

Having given some thought to this, I'd still use machine learning. I would attempt to augment the negative controls however. Is there a suitable published data set? R program is caret but it high barrier to entry without appropriate support.
250 samples is just about enough. If you can augment that to 350 with negatives ... should be enough. What you do is dump all your data ("features") into a single record, including the methylation data and use toxic element exposure as the training target. The output you want is called "feature selection", what you are wanting is methylation sites to carry more "weight" than any other feature. Then you start other analytics to see what it's associating with. Its not trivial getting it working and a sceptic would complain about the lack of significance in the results, but its a good starting point.
GWAS might work, just considering methylation vs toxic element exposure. SEM, I've never used it so I don't know.


text:
Exporting Panther functional classification of gene list
After inputting a list of gene and running an analysis in https://www.pantherdb.org/ (such as "functional classification"), I get a result table that is limited  by an X number of entries per page. I can select some entries from the table and export that subset as a text. But, Is there a simple way of exporting/downloading the whole functional classification for my gene list as a txt?

Do you mean this? https://i.imgur.com/ZF7C06f.png Or is that what you tried and it doesn't work for you?

I tried that but only exports whatever is paginated on the html OR whatever I can select on the page. I want to get the table for all the thousands of IDs I uploaded.


text:
Estimating rate of evolution from a phylogeny
Any suggestions for the best way to estimate a metric for evolutionary rate of viral sequences of the same species?
I was using terminal branch length from a phylogeny and would quite like to try sub rates but i am really new to this

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.