text:
How do I calculate DEG cutoffs?
So for context, I have a set of TPM values for multiple genes for different samples, and I need to calculate the differential expression for RNA-seq values. I've been using this website as a guide: https://raw.githubusercontent.com/ucdavis-bioinformatics-training/2019_March_UCSF_mRNAseq_Workshop/master/differential_expression/DE_Analysis.Rmd. The overall aim of this calculation is so that I can use PANDA to map TM motifs/combine gene expression data to create networks representing interactions between transcription factors and genes, and GSEA to analyze genes ranked by fold change or differential expression p-value.
Note that I am using log2(TPM+1) values. How do I calculate what the filtering cutoff should be? Note that I already excluded all genes that had any values of 0 for any of the samples, and I will use the Benjamini-Hochberg method to adjust p-values for multiple testing afterwards.

The number of observations is little less than 19,000

TPM + Limma looks fine.  Cut-offs are cut-offs, end of story.  The DEG issue is the variance is much greater than the mean: DESeq2, edgeR and voom attempted to "hardwire this relationship" hence criticisms from Li et al (2022) and others. DESeq2 and edgeR need tests to assess the fit between data and model (residual). That's common in other stats.

TPM per se has been criticised in some published comparative studies  recently. I wouldn't discard your approach because: 1) it's a recognised approach; 2) transformations are complicated and the alternative (counts) are back to edgeR, DESeq2, NOIseq ... not fun; 3) the weaknesses are not formalised like 20+% FDR (theoretically always anticipated for DESeq2 etc.). You may be requested to do an alternative method post peer-review, or you could go into comparative methodology (which is good analytics). Your approach is ok stick with it, complimenting it is your call

tl:dr use limma (not voom, that's for raw counts) and refer to previous questions:

https://support.bioconductor.org/p/98820/#98875
https://support.bioconductor.org/p/125144/#125156

The problem you have is not so much the threshold but the method to infer any metric to be used for this threshold. The arguably most common methods in R for DEG analysis are DESeq2, edgeR and limma, based on number of citations. The former two require raw counts which you apparently do not have, but if you do then use this with any of the mentioned tools.
limma is what is left and could be used here with TPM. Note that TPM has the problem that counts are divided by a length factor. A brief background: Any statistical test aims to decide whether the given variance between two groups (basically the logFC) is greater than expected by chance, which is the variance within group. In RNA-seq the variance is not equal across the range of counts (low, middle, high, for simplicity). On log2 count scale variance is greater for genes with low(er) counts than for higher counts. Usually, using raw counts the method (DESeq2, edgeR, limma, and others) learn from the data what is the expected variance at the given average expression level, and incorporate that into the model to make adequate inference even at small sample size. That is by the way where the common assumption "most genes do not change" comes from. Assuming that most genes are not DE, then the mean-variance trend represents the expected variance at the given average count range, and this can then be used to help decide whether a gene is DE or not. With TPM you basically throw this out of the window because the actual counts (as they come from the sequencer) are manipulated by dividing by length and the entire mean-variance model gets distorted. Another factor is that TPM, which corrects by total library size can fail to properly normalize counts, I wrote a summary here:

https://www.biostars.org/p/9465851/#9465854

Long story short: TPM is not good for DEG.
If it's all you have then use what is called "limma-trend", see the limma manual at Bioconductor, which will directly model the log2-TPMs and give some logFC and FDR. You can then use standard cutoffs, such as FDR < 0.05 & logFC >/< 1.5. These cutoffs are not perfect, but common, and given that TPMs are bad anyway you can and should stick to standard cutoffs as much as possible, rather than infering your own ones. This at least gives some standardization to this suboptimal analysis.

Two solutions - the underlying issue is called Type 1 error (below).

If you are performing studies on humans then the Benjamini-Hochberg method is without question the method of choice, because this includes the false detection rate for EdgeR and DeSeq2. Thats been independently calculated.

Thus, the paper that caused a stir was Li et al (2022) who wrote,

Particularly, the actual FDRs of DESeq2 and edgeR sometimes exceed 20% when the target FDR is 5%.

Benjamini-Hochberg (BH) method solves that for humans (EdgeR and DeSeq2), because the problem with the BH method is you need a false discovery rate, Li et al 2022, calculated that its about 20% ... so for EdgeR and DeSeq2 the problem is solved. The BH tests will simply use the estimate for the FDR by Li et al (2022) instead of the FDR by EdgeR and DeSeq2. Cool? It is :-) it wasn't what Li et al (2022) intended but it is the solution to a difficult conundrum (concerning a statistical issue or parametrics vs. non-parametrics). So yeah it's cool (needs count data).
To perform BH you simply rank the observed probabilities (e.g. from Limma), and there are significant if each is below (or equal to) the BH rank. The BH rank is the total number of tests as denominator, the rank position as numerator and this is multiplied by the false detection rate as a frequency.
HOWEVER, you are using Limma and I don't how the model works. I know the underlying stats behind EdgeR, DeSeq2 and NOISeq.

Okay so the false discovery rate is not known ... in the way Li et al (2022) calculated it. What now? Limma is likely to calculate an FDR, you could use that via BH. My system is so far from humans there's no way Li et al (2022) FDR is relevant.

Bonferoni's approach - I don't know how mainstream this is, but its what I use (I've a background in classical stats), because I've no idea what the true FDR is. This approach can be exemplified via Volcano plots is to visualise the Type 1 error issues.

Here's a volcano plot which permits online viewing via links obtained here. I can't change the size BTW (Bioconductor and loads more BTW don't permit this),

This is using a flat probability - which in the picture is very low - it's the horizontal threshold (y-axis). So the P value is one parameter approach to significance ....
The other parameter is the log2 threshold ...
Here's a two-parameter approach obtained from GeneKitR's site, who do permit images to be link viewed (thank you genekitr).

Here the authors have used log2 as the threshold (x-axis). The threshold is probably log2 > 1 and log2 < -1, which is reasonable. They have set the P-value threshold to be really low, looks to be about 5%. The plot reveals the effective P-value is 1E-2.5, which is low, but nothing really starts happening until 1E-3, which ... okay. Although this is a two-parameter approach the guys really just relied on log2.
Type 1 error and a 1-parameter method threshold
The massive problem is the shear number of tests performed and this will result in significant values PURELY BY CHANCE ...artefacts of critical probability. BH is one approach, Bonferoni's is an alternative way in.
Put simply for your data to calculate the Bonferoni threshold if a 5% critical probability is sought then its ...

0.05 / 19000 = 2.6e-6

Which is okay. So any single result <2.6E-6 is significant. This okay, but you will lose results so its vulnerable to Type 2 error. Why do you lose stuff ...
Type 2 error and 1-parameter method threshold
The biological question is e.g. is 0.1log2 biological sensible? What does it actually mean? Honestly I don't have a clue ... and the problem with a 1-parameter approach, it's very conservative but infers weird stuff. An expression difference of 10% isn't really serious, but thats what it's measuring. ... because its measuring this it means the probability via Bonferoni's is set high - for meaningless stuff and you can be missing important shifts in log2, e.g. >1 because they are falling under the high threshold.
Type 1 error and a 2-parameter method thresholds
This is where it all starts getting interesting. As soon as a +/- log2 threshold is set the number of total tests will drop A LOT. That means the P-value threshold raises. For example, if selecting > (or <) +/- 1log2 as used in the above example, supposing the number of observations above/below the +/- 1log2 was 5000. The P-value now raises via Bonferoni, e.g.

0.05/5000 = 1e-05

That looks better there'll be more results. Supposing we set a threshold of  +/- 1.5log2 and obtained 500 results ... then the new Bonferoni derived threshold is,

0.05/500 = 1e-04 (0.0001)

That looks fabulous giving loads of significant results ... Wait isn't this all getting a bit crazy, because it depends on the biological assumption of what a reasonable log2 is. Yep, indeed, agreed.
Type 1 error and a 2-parameter method thresholds: weaknesses
The issue is that a reasonable log2 is sensible, because there's a point where log2 expression differences becoming utterly meaningless. However there's a point up to -/+ 2.5 log2 (which is very high) where you could choose anything - that shifts the critical probability and that shifts how many results are positive.
So how do you ...
Set a log2 within a 2-parameter threshold?
I've just solved that. I was really happy. The practical answer is it's complicated. My contractual basis is I can't reveal that because public disclosure isn't permitted without top-level approval.
The approach is to use an independent method to calculation DEG  (if you check my posts you'll see I do AI). You can then perform parameter tuning of DEG against the independent method and statistically assess the association. In my case I was surprised by the answer - intuitively I wouldn't have gone there, but my system ain't humans.
Conclusion
In your case,

simply check the literature to see what other log2  thresholds are in use, use that and adjust the P-value threshold accordingly as described above and cite that work.
look at the biological reality of the results at different log2 thresholds. For example, an operon should behave like an operon. If the log2 threshold is smashing-up that operon - that doesn't make sense and adjust log2 threshold by lowering.
alternatively use BH, this however assumes Limma's FDR calculation is correct and the problem is there no independent study examining that.

Either of the above 3 approaches will be fine your context. This is not cutting-edge statistical methodology it simply needs to be reasonable.
Ultra-conservative parameters and Type 2 error
To avoid all the complexities some investigators simply use > (or <) +/- 2.5log2 and P-value = 1e-6. These values in ANY system are extreme when used together, but there is no chance of Type 1 error. The danger is Type 2 error , where significant values are either omitted or considered non-significant. This could happen here from two fronts:

If this results e.g. 100 observations, and DEGs with P-values < 5e-4 but > 1e-6 are omitted
DEGs < (or >) +/- 2.5 log2 and above the lower threshold are omitted.

However, it will depend on the question.
Note Much of the theory in point 2, would also apply to the BH-method. However, I've worked this theory exclusively for Bonferoni's correction and I've not given much thought to BH. This is because it's not relevant to my purpose.

See comments below your post about your TPM approach. This post is cut-off criteria.

Just to finish this post, here's a volcano with thanks to Britannica  for permitting the link-view. I think the idea of a volcano plot refers to the "plume" of ash, dust and occasionally lava (which some of us have real-world experience).


text:
ZIP the fasterq-dump files in loop
I want to zip the files generated by fasterq-dump simultaneously as they are generated. The code that I tried is this.
cat SRR_Acc_List.txt | while read i
do  fasterq-dump ${i} --split-spot -v --skip-technical  -e 15| gzip -k  ${i}> ${i}.gz
done

Haven't tried, but I believe you may add -Z to fasterq-dump and remove -k ${i} from gzip. In addition, I would recommend the following:
cat SRR_Acc_List.txt | xargs -i echo fasterq-dump {} --split-spot -Z -v --skip-technical -e 15 \| gzip > {}.gz | sh

It is simpler and more flexible. For example, replacing sh with parallel would allow you to download in parallel.

I tried this and it worked
cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-spot -v --skip-technical  -e 20 ; for j in i ; do   gzip ${i}.fastq  ; done; done

source: Stack Exchange

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

sort: votes


text:
How to handle technical replicates in cellranger count?
I'm working with cellranger count to annotate some publicly available fastq files (SRX7888073). For each sample, 4 runs are available, which I assume represent technical replicates.
How do you recommend to handle them in cellranger? Are they 4 different samples (i.e. 4 separate runs), or 1 sample with 4 lanes?
Thank you !

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.

Welcome @David, nice you joined us. It could be helpful if you could point to the information on the site, ie. "Male 3; Schistosoma mansoni; RNA-Seq". Thats kinda key, clicking offsite links to answer questions is sort of discouraged. The actual experimental details will be explained on the site ...

It depends what you mean by "technical replicates".  If it is the exact same library prep spread out over two lanes, or two flowcells, they must be combined at the fastq stage.  Cellranger will naturally do this with fastqs that keep their original Illumina-assigned names, but just have the files split by lane.
If one file has a read from a gene with a particular cell barcode and UMI, and another read in a separate file but the same library prep has a read from the same gene, with the same barcode and same UMI, those two must not be both counted. So if you combine the counts, that's wrong.
A NextSeq naturally has 4 lanes, so that's almost certainly what you are looking at.  The read names would confirm that, unless they've been altered.  Cat them all together, or alter their names so that they match the original Illumina specs with lanes 1-4.  You might have to do that anyway, cellranger used to be very picky about that, and might still be.

Agreed on all of that. If you alter names though then furst make a new folder, symlink files into that and alter the symlink names. Don't  alter names of actual files, that's error-prone.

If they are truly technical replicates in the most specific sense (e.g. from the same extraction and library prep, just split up for sequencing), then the counts for each identified cell can be combined (i.e. added together). One way to do this is to concatenate the files prior to running any data processing scripts.
If there is anything different about the runs, then they should be treated as non-technical replicates and properly integrated:
https://satijalab.org/seurat/articles/integration_rpca.html
We had a situation with Rhapsody single cell data where a library was split into three on three different cell barcode plates, and was showing a small amount of library preparation bias that influenced clustering. Carrying out an integration of the libraries helped to reduce this bias.

Catting the fastqs will not be the same as combining the counts.  If the catted fastqs are processed by cellranger, deduplication will be handled properly.  If the fastqs are processed separately, you will overcount, because duplicates between runs will not be accounted for.

After some digging, it seems like each of the 4 runs has I1, R1 and R2 files - they're replicates.
You can either

run them separately and then perform cellranger aggr; or
link to all 12 files using soft-links that follow a proper naming convention (the files need to have the same prefix, lane numbers need to be different per replicate), ensure all softlinks are in the same directory, then use the directory with cellranger count.

I'd recommend the former approach as it will probably save a little compute time, is parallelize-able and will be easier to debug/drill down into.

using cell aggr is totally wrong  If you had 4 different samples, it would be a coincidence if some of them happened to have the same barcode, and the right thing to do would be to treat those as two separate cells.  With this, if the same cell barcode turns up in each lane, it's actually the same cell, and they should be combined.

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


text:
Variant Filtering Using AnthonyOzerov/variant-filtering Github Repository And Python Errors
I am using the following Github repository that I cloned. One apparent issue is that it is using an old/defunct link for the phenotype_to_genes.txt that is no longer active and after several tries I found what seems to be the correct link. However, I cannot seem to get it to work as it seems to be throwing up an error that makes no sense to me when I try and run it via a bash script:
Variable I defined in advance:

pedigree_file has the path to the pedigree file with 4 columns (Family ID, individual_ID,Status [Child,Father,Mother,Sibling], Sex [Male or Female] and Phenotype [Affected or Unaffected]
Mendelian a path to a file created via Mendelian filtering
output_dir path to the output directory
phenotype_path a path to a file with 2 columns group/family and then HPO
VC_Python_directory path to the variant filtering directory
Map_File path to the phenotype_to_genes.txt
family_names a list of family/group names

The bash script:
cd "$VC_Python_directory"

IFS=',' read -ra family_array <<< "$family_names"
for family_name in "${family_array[@]}"; do

    # Concatenate the file path
    file_path="$output_dir/$family_name.txt"

    # Echo the family name and file path
    echo "Family Name: ${family_name}"
echo "File Path: ${file_path}"

    # Run the Python script with the correct arguments
    python main.py -p "$pedigree_file" --data "$Mendelian" --output "$file_path" --family "$family_name" -ph "$phenotype_path" --mapfile "$Map_File"
done

The error:
Family Name: Sample_1
File Path: /path/to/output_directory/Sample_1.txt
Path/to/variant-filtering/utils.py:77: DtypeWarning: Columns (2) have mixed types. Specify dtype option on import or set low_memory=False.
  phen_to_genes = pd.read_csv(mapfile, sep = '\t', header = None, comment = '#')
Getting relevant genes for family phenotypes...
Traceback (most recent call last):
  File “path/to/variant-filtering/main.py", line 26, in <module>
    load_phen(families, args.phenfile, args.mapfile)
  File “path/to/variant-filtering/utils.py", line 80, in load_phen
    phen_to_genes.columns = ["HPO-id", "HPO label", "gene-id",
  File “path/to/python3.10/site-packages/pandas/core/generic.py", line 6002, in __setattr__
    return object.__setattr__(self, name, value)
  File "pandas/_libs/properties.pyx", line 69, in pandas._libs.properties.AxisProperty.__set__
  File “/path/to/generic.py", line 730, in _set_axis
    self._mgr.set_axis(axis, labels)
  File “/path/to/managers.py", line 225, in set_axis
    self._validate_set_axis(axis, new_labels)
  File “/path/to/base.py", line 70, in _validate_set_axis
    raise ValueError(
ValueError: Length mismatch: Expected axis has 5 elements, new values have 7 elements

I have checked over the pedigree file and map file and phenotype path files but I am really unsure of how to go about this pipeline as I have never seen this kind of error before. I am not sure I can ask this here but does anyone perchance have an inkling as to what is causing the rather peculiar error?

Apparently the format of the phenotype_to_genes.txt File changed, while at the time when this pipeline was published the file consisted of 7 columns, it only has 5 now. That's what the error message tells you. It's best to file an issue at the github repository to get in touch with the authors. 1) so that other potential users are also aware of this issue and 2) the authors might be willing to fix this issue/help you get an updated version of the pipeline or an older version of the phenotype-file.

I was wondering what the original has as I might be able to publish the changes but I doubt I can get in touch with the author again. I think from what I can tell the original is 1. HPO-id: 2. HPO label:  3. gene-id: 4. gene-symbol:  5. gene-name:  6. gene-source:  7. disease-source:

Thank you @Patrick.


text:
Converting convert .svs files into .mrxs format
I wish to convert Aperio AFI microscope .svs files from a pathology lab microscope into 3DHistech .mrxs format. Lab ppl have some pipeline X and they asked me to look into the matter.
I can read the SVS files via openslide or openmicroscopy. Conversion to 'regular' pictures (tiff, for instance). I figure out that since openslide can read mrxs, writing this format should be feasible too, but can't find any (free) software to do this. I don't care about any specialties, just need the end result to be readable by Histech readers.
Can anyone help with the conversion? Lacking full software, I could try to code writer if resulting file description is available.

Can you please add more context around your actual problem? I'm particularly interested in why writing as .mrxs format is necessary. What do you want to do with it after the conversion?

It's just that it's unusual for conversion to a non-free format to be necessary, and can be quite a difficult process to get it right. Understanding more about the problem may help to find a better way to go about doing it.

I've come across this problem before - Aperio SVS files are difficult to wrangle, and are best dealt with using proprietary software.

I had a look at slideio (python module) and this will happily import SVS as an image object, so e.g. matplotlib will work on it .. cool. However there is very little supports for MRXS. I know this because of expert articles such as this, https://www.openmicroscopy.org/2016/01/06/format-support.html which is part of openslide. There are saying they wouldn't do it. Certainly the question isn't solvable on this site. Summary, SVS is easy BUT not exporting to MRXS. As @RamRS points out some sort of TissueGnostics purchase looks the only option. Commercial companies keep you locked-in. Upvoted.

If you truly "just need the end result to be readable by Histech readers" the simplest (and maybe even cheapest in the long run) solution might be a paid ($) license for CaseViewer (or whatever their newer viewer is called). I have no idea what that would cost but many of these vendor wsi viewers can read and even write other formats for a price.


text:
Pacbio HIFI pbmm2 alignment metrics
I am new to pacbio sequencing data. I just did some alignment of pacbio HiFI data using pbmm2. I have the bam now, I would like to collect some metrics on the alignment. I used to work on the illumina data, I used picard CollectAlignmentSummaryMetrics, CollectGcBiasMetrics, CollectInsertSizeMetrics to collect a series of alignment metrics.
I was wondering if there is similar tools for pacbio data. I searched and did not find an answer.
We are mainly interested in aligned rate, sequence length, coverage, accuracy, and GC content. So I write some scripts for it. These scripts are inside the nextflow, the input is the bam file
samtools flagstat $bam \
 | head -8 \
 | tr " " ","  \
 | tr "+" " " \
 | awk -F',' '{ print $1, $3}' \
 | paste - - - - - - - -  \
 | awk '{print $1, $3, $5, $7, $9, $11, $2, $4, $6, $8, $10, $12 }' >  align_prp

 samtools depth $bam \
 |  awk '{sum+=$3; sumsq+=$3*$3} END { print sum/NR, sqrt(sumsq/NR - (sum/NR)**2)}' > coverage

 samtools view $bam \
 | awk '{ $35=length($10); sum+=$35; sumsq+=$35*$35} END { print sum/NR, sqrt(sumsq/NR - (sum/NR)**2)}' > length

 samtools view  $bam \
 | cut -f 10 \
 | fold -w 1 \
 | awk '($1=="G" || $1=="C") {N++;} END {print (N/(1.0*NR)*100);}' > GC_content

But I am still interested to know if there are some package available there.

Is this useful at all? https://www.biostars.org/p/9494195/

No, I am not accessible to Revio, Sequel II and IIe systems.


text:
How do I fix this limma line?
I'm trying to conduct lmFit on my fread::table matrix:
fit <- lmFit(matchedgeneTPM)[,-c(1,2)]

It's giving me the error that expression object should be numeric and that there are 2 non-numeric columns (because the first two columns are the labels of the genes and gene ID's I'm working with). I thought the last portion of my code would allow it to skip over the first two columns but apparently it's not the right way to do it. Can someone show me how to format the code correctly, please?

My guess is you need to index matchedgeneTPM before you run lmFit(), not after:
fit <- lmFit(matchedgeneTPM[,-c(1:2)])

That worked, thanks!


text:
how to search and calculate occurence of certain word (like TAA) in fastafile
>003197.2:c1179912-1178219_reverse_complement
AGGCTTTTAAAAGCCTCCTGGGTTTTTAGTGAAGCTGAGTGATAGAAGCTCTGTATTTGCATAGCGTTTT
GTCGGGCGCTTTAGCCGGCGCTTTGCCCTGGCCTGAGAGAATCTGCATTCCGTTGTATAAGGTTTTTTGT
CCTTCTGATGCTGTAGATAATTCCCCCATGTCGCGCCGGGTTCTCGCAGGACAATAATTTCCGGGCGAGC
AGAACGGTTGCCGCGACGGTAAGAAGATCGCGCTGTAAGTTATAGAGGTTATGCAGCGAGTGGTTAGACG
CCGGTCCGCTTTAACTTTGGCTAACTCCATTTGGTTTGCCATTGACGTTAGAACCGGGTCTTGTTTACCC
GTGCTGCAATAAGTTCGATAAGATTTTTCTTCAATGCTTTTGCCGCGGCTTCTTCTTGTTTTGTTGCTGG
GCATCTCTAAACGCTACTGCCGCAAAGCGATGAGCTTCTTTTGCAGGTAAGCCATCCTGCTGCTGAGTGC
GTGATGCCCGTTATGCGTGAGTGTATTTTTTATGGTTTGCCAGGGCTGGTTATTAAGCTGCTTGACCTGA
CATAAGCACTGGGAAAGATATCTTTTGCGCCGATTTTCATCTCTGCTGCAGGGAGCTGCGTGTTGGTATA
ACCGTGGACATCCACAAATTATTGGCGTGATGAATATTCTTGGTATCCCAACTGCATACGCCCTTTCCCT
ATAGGGGGAAAGCACGCCATGACGTATCCCGCAAAAAAGCGTTTTATCTTTACCGTCCTCATGCACACTC
GTGCCGCAGTTAATACTTCTTTGGCTTTGTTTTCAGCGCCGACGTGACGCAGAAGCGGATCTTTTTCATG
ACGGATACCAGTTTCAGGCTTACCGCCTCGCCCGCTAAGGCTTTGTTAAGCAACTCAGGTTTACTAAAAA
CCATGCGCGCATTTGGTCCTCGACCATCGTTCCCTCTTTGCCGAAAATATTCGACGCGGTGAGTAACCCG
CCGTCTGTAGATCGCCATCTTTATTGCGGATTTTTAAATGAATCATTTTTCCCGGCTGGGTCAACGATTG
CCAAAGCCGAGCTTGAGCGCCAGCTCATTAACACCCACATTAAATGCGGCGACGTCCGGTTTTATTTTTA
TTCAGGGCGTAAATCATTGCCTAATAACTGATGTAGCGCCTCGGCATTATAGCTATCCGATGCCTTAAGG
ATGTATTGACGACCTCATAATTATCCGGGTATTGCGCCAGCCATTCGCCAACCCAGCCACCTGGTCTGGC
AGTTTATAGGGTTCGCCGCCATCTTTATGATGTTGGTTATTTTTCCATATATCCTTAATCTGGCGCGCTA
TTTGCCGCTTTTACAATTCCAGGCGGGTACCGCGTCAATTTCATGGGCTAACATGGCAAGGCGTTGTGCG
TTAACATATGGGTCTGATGTAAGGAAATGATCTCTCGCTTGATTTCTGAATCCATCATCCCTGTACGATC
TTACCGCTATTCAGTAATACTTTTTGGAAAATTTTCTGTCCACCGCTATCCGGAAGACTACCAGGCGCAC
CACCTCTGGCGATAAATTTTTCATTACTTTGTTTCCCGCCCCGCCCGTATTTTGTTTCTGAATCTCCAGG
ATGAAGAAATGCCTTTTACTGACTGCCAAATATTTTCATCCCCAACTCGTTTTTGATAGGAAAGATTGAG
TCAAGATGTGATTA
>003197.2:c3026750-3024686_reverse_complement
TTTGATCTCGGTCTGCATACCTGGCATTATGACGGGGGGCTGAGTCCTTACACTTGTAACCATTATTAAT
CCCCTGACAGCGTCGCTGTGGCACTTTCTCTGACTGCGGAAAGATTCGCCGCAAGATTCGTGGCCTGCGT
TTAAACAATGCATCCAGACTCGCCTGTTTGATCAGCGCGGGAAAATCTTCCAGTTGCGGGCCTTTAATTT
CGCTTTCTTACCGGCGACATTATTTGAATTGGTAAAAACTTCTTTTAACGCCTCAGCGTCTTTCCCGCAT
GCGTCTTCGCCTCAGGAGAATCACTGGTGGCGTTGAGCGCTGAACGAAAGAGCCCGGCAAACTCCATTAT
AAAGCTGACTTTTCTTTCAGGCCATCTTTTATGATTTGCGCGGTATACTCTGCCCCCACCTTCATCAGTA
ATCTAACAGCTGCTTTTCCAGGTTTTCCAGCTTCGCTTCCGCTTTCTTTGTTTCTGGCAGCCATGGCCCA
GCATCACCTCCTGCATAACAAGATTCGTCGATAATTTACTTAATTCACCGCCAGTGTTATTTTTGATAAT
GTAAGATTGCTGCGGGTTAACGGATCTAACGTACAGCCAAAGTTATGTTCAATGCAGCTGGCAATATAGG
TTGCGTTAATTTTTGGTGGCACATATCAAGCGCCTCAACAGCCTTCGCCGCCGCTTTGTCAACAAGGTGC
GTAGCAGCAGGGGTATCAACGTTTGCATTTCAAGGTGCCGGGCTTCCCGTCCTACGCTGGTACCCTGCTC
GGTCCAGCCGCAGGCTCCGGAATTTTAGGGTCAGGCAGTTTATCTGCAGGAATTTGTGCAAAAACATTAC
GTTACTGCTATCAATATTAATATTGATATTAATACCCGTAGGTTCTGCTTTTTTCCCACCATCAGGTACT
CGATATGGCGCTGGCTATTATCGACATGGCTTCGGCTATTGTTAATATGCTTACTGTTATCCACGCTATG
CGTTGGCTGTTATCAATAAATGTTCGGCTATTATCAATCGTCTTCCGGCTATTATCATGGTTGCTATTAT
AGTCGAATCCACACGCGAATGACTATGACTTACGTTGCTGCTGTTAGTCGAATGGTGTGACTCGCCATTG
CAATTTTACCCGCGATGCCATGGTCAATTGCCCCTGTGCTGGCGCTGTGTGCCGTCTCCGTTTGATGCGT
ACCACCTTTCCATCTTTACTTTCGGATGAAGCGTTGGTCACAGCCTCTGCCGTCGCATGAGCCGTGACGT
GCTCTTGCTGGTTACGCCATCGACTTCATCAAAAGAGGTTGTTTCACCCGTAGTGCCCTTTTCTGACGTG
CTGCCGTCTGGCTTTGCTGTTTATTATCATCAACGCCGTGCACCGTGGCCTGTACCGGCTTGCCGATAAT
TGTAAGGTATCGGTTTTAACATTCTCGACGCCAGCCAGTTGACTGGCAAGCGGCTTCACATTCACAATCT
CTTATCTGTCTGACTATTATCATTATCGGTCGTACCGGCTTTATTACCTGTAATTACCGTCGTGTCTGAT
TATTCGTAACGCTCGGTACTGTTGAGAGGAAACTATTTTGCTTAAACTTCAAACCAGAAAATGGGCCAGT
GTATCTGGCTCCAGGGCTTTACGTATCACACCCAAAAACGTCTCACGGGCGTCGAAATGCATTGAATGCA
TATCGTGTCATTGGGTAATATTTCCGCACGCAAGCCGTCAAACGCTCTGCGTACAGGGAACGGTGTGGAG
TCTCCATTTTCGCTTTCAATTCTGGGTGCTTGTCAATATCGCTACATTGCGCTTTCAACGCTGCACTCTT
CTGGCAAACTGCAGTGCAATCTCAGCCAGTTTTTCTTTTTGTGGATGATGAGTAATGACCTCTTTGAGAG
GATAAGGCCGTCCAGCACATTACTTAACACATAGTCAGTTTCCCCTTTTAGCCTGGTCAGCCCCGCCTCT
TAACCCGATCCACGCCAGGTTTATTCAGGTATGACTCGTAAGCCGGGCCAGCTCGCCAGCTACCGTCTCC
TTAACGCTGCATGTGCAAGCCATCAACGGTAGTAA
>003197.2:c3010816-3009905_reverse_complement
TATTACGAGTCGTAGGACTTAGCATACTTTTCCCTCCACATGATAGCTCCTGCACCGAAAATATCATCTT
ACAATACAGGTCAATGCTTCGGTGTTTAGCTCCCCGCTGACGTCAGGCTGAGGCGATAAACTCGCCCCCA
TTCGAGATCGGTTTCCTCGTAACTGATATGAATCCAGCTGCCATCTATAATTTCACTTTCCAGACGCTCA
CTGACATAACAAATTTAAGATTTAACTCTGGATACTTCTTATTGGCCTGTGCAACAAGAAAAGGCATCAT
TTTACGAGAAATCGGGCAGATTCAACGCCTTCCATTCTGGTCTCCTTTATAGAGGAAACAAGCTCATGGA
AATCACCGATGTCTTTCCGTCCATAACTCTAAAATCGACGACTGAAATATGGATACCTGAACTTCCCATA
CAAGAGCTGCTTTGGTCCTCAACGCCAGTAAAGCAGGTCCAAAAGCGCTACACGCTGCTGGTTCGAACAA
CATTCAGAAGAGCTTCGTTGAATGTCCAGCTCGACCATAGCAAAATAACAATCAGGCAGTTGTTCACGTT
TTCATGAATTTTTACTAAGTTCATAAATTCAAGCTGAAGTTTTTTGGCGAGCGCCAGGCTAAAAATACCG
GCAGATAGCGGTCTGCTTTATCGGACGGGAGAAAAGGTTCTTCACCACACAGACGTTCACAAATATTATC
GCCGCCGGATTGGCCTCCACATATTCATTTAATCGTTGTGCGCCTTGAGTATGTTTGTAAAAACTCACCG
TAGCATAACGGCATTGTTATCGAATCGCTCATAAAGCGTTTCATTCTTTTTGTTTACTATACTGCTTCCC
ACGGTTTAAGTAAAGACTTATATTCAGCTATCCTTTTTTTATGAGCGGATATAGAGAGTTTTTTATCGTT
TT
>003277.2:c27903-27178_reverse_complement
CCCCATCATAAGCAGCTACAATATTCAAAGGAGGGATGTTTAGATTTAGATTAGGCCTATTTATGGGCAT
ATCTTCCGCATTTGGTTGTGCAAGGAGTTGAAATTATTTTTCAGGTGCTTATTTGTAGATGGTATTTCCG
AAAGAAACCACTCTGGCGCATCCCTGAATAGTCAGGCACATCCAGCGCCTCTTTAAAGTGGGATACCGGC
GTGATTCTCGACGAGCATGGATGAAAACATCGTAACCATGGTTATTCAGCTGAAAACCTTGGCTCATAGC
GCTTGTGGCACCATATCCCTGAGCACACTGATGTGGAACTTGTCACCGGCAAATTTGCCCTGAGACTGAG
CATATCGGTCACTTTCCACTTATCTACCGGACTGTCCTCTGAAAACAGCAATCCGGACAGCGCTTGAAAT
CCTGGTCTGGTTTTATATACAACGTGAACTGAGCGCCCAGGCTAACACGGGCTTGTTGAACGACCTTCTC
AGTCTGGATTCCAGACACTCTATAAATTGCCGTGTCTTGTGGAGAAACGACGCACTGTACTGCGAATTTT
GAGATATTTCCAATTTTCAGGATGAACGTCTGACTCAGGACACTGTCCTGAAATAACCCCATTTTCGGAT
GTTCCTCACGTAAAGCCTGTCTCTGCATTTCGCCACCATCACGCCCACTTCGTAGTTCATTACGATAACT
TTACTCTGTCATCAAACGATAAAACG

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.

Welcome to the site @tonu. Please make your question clear within the question. Finally, what coding are you happy with, bash(shell), R, Python, Perl, awk etc ...?

@Tonu this has been correctly answered could you please consider upvoting/accepting?

Since you tag python, I suggest a simple Biopython solution:
from Bio import SeqIO

print(sum(s.seq.count("TAA") for s in SeqIO.parse("file.fa", "fasta")))
# 97 on your input file

If you wanted an overlapping sequence count, Biopython has a method for this too:
print(sum(s.seq.count_overlap("TAA") for s in SeqIO.parse("file.fa", "fasta")))

Here they are the same, but for a string like "TAT" they could be different:
>>> from Bio.Seq import Seq
>>> Seq("TATAT").count("TAT")
1
>>> Seq("TATAT").count_overlap("TAT")
2


text:
How to find the correlation between two genes and their separate overall survival outcome?
i am struggling to find the correlation between two isoforms and their survival outcome. i do have a dataset that is a cohort for thousands of genes and am only interesting in two which they should be isofroms to each other. However, i thought the ratio between these two might be a prognostic marker. i need to find a way where i can see if this is possible by interpreting the clinical outcome which is the overall survival. many thanks guys.
Crossposted: https://www.biostars.org/p/9566172/

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.


text:
adding multiple files in folder and zipping in loop
Through fasterq-dump, I am generating 2 files in 1 run of loop. I want to club these two files (A_1.fastq & A_2.fastq) in one folder and zip the it.
cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-files -v --skip-technical  -e 15 ; for j in 2 ; do   gzip ${i}_j.fastq  ; done; done
cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-files -v --skip-technical  -e 15 ; for(j = 1;j<2;j++) ; do   gzip ${i}_j.fastq  ; done; done

cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-files -v --skip-technical  -e 15 ; for j in 1 2 ; do   gzip ${i}_${j}.fastq  ; done; done

cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-files -v --skip-technical  -e 15 ; for j in 2 ; do   gzip ${i}_${j}.fastq  ; done; done
cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-files -v --skip-technical  -e 15 | mkdir {i} | for j in i; do   gzip ${i}  ; done; done

cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-files -v --skip-technical  -e 15 | mkdir ${i} ; for j in i; do   gzip ${i}  ; done; done
cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-files -v --skip-technical  -e 15 ; mkdir ${i} ; mv ${i}_*.fastq "${i}"; for j in i ; do gzip ${i} ; done;done

cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-files -v --skip-technical  -e 15 ; mkdir ${i} ; mv ${i}_*.fastq "${i}";done; for j in i ; do zip ${i} ; done

cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-files -v --skip-technical  -e 15 ; mkdir ${i} ; mv ${i}_*.fastq "${i}" ;  zip ${i} ; done

cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-files -v --skip-technical  -e 15 ; mkdir ${i} ; mv ${i}_*.fastq "${i}" ; tar -cv ${i} | gzip > ${i}.tar.gz ; done

None of them is working.

Try your code on one set of FQs and use a loop only after that set gets processed successfully. Also, why pipe tar to gzip when you can simply do tar -czf archive.tar.gz files_to_archive?


text:
Why is my centred log ratio transformation on a distance matrix in R not symmetrical?
I am trying to transform a distance matrix from the qiime2 pathway using a centred log ratio in R. I uploaded the matrix directly from a qiime2 output file, and have converted to a data matrix and have checked that it is symmetrical.
structure(list(X1 = c(0, 0.2177609998, 0.2133361674, 0.1549136105, 0.1400395799), X11 = c(0.2177609998, 0, 0.07805820645, 0.1418994689, 0.1934668819), X12 = c(0.2133361674, 0.07805820645, 0, 0.1475390242, 0.1857477705), X13 = c(0.1549136105, 0.1418994689, 0.1475390242, 0, 0.1046740994), X14 = c(0.1400395799, 0.1934668819, 0.1857477705, 0.1046740994, 0)), row.names = c("X1", "X11", "X12", "X13", "X14"), class = "data.frame")
dm <- read.csv(file = "dm.csv", header = TRUE, row.names = 1)
isSymmetric(as.matrix(dm)) [1] TRUE
dm_matrix <- (data.matrix(dm, rownames.force = NA))

I then tried to transform the dataset using the dm_matrix_clr <- data.frame(clr(dm_matrix)) function from the 'compositions' package. This produced an output, however it was not symmetrical. From my understanding, the transformed distance matrix should be symmetric like the input matrix. The same numbers on opposing sides of the diagonal are being transformed, but are resulting in different outputs.

Any help on how resolve this problem so I have a symmetrical transformed distance matrix would be greatly appreciated.

as i understand it, the transformation assumes compositional data (columns sum to 1). Distances are not compositional. Where did you get the idea for this transformation? Can you link anything? Not really familiar with qiime or surrounding workflows.

Here's a paper with a good explanation: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6755255/
"All components in a composition are mutually dependent features that cannot be understood in isolation. Therefore, any analysis of individual components is done with respect to a reference. This reference transforms each sample into an unbounded space where any statistical method can be used. The centered log-ratio (clr) transformation uses the geometric mean of the sample vector as the reference"

that explains why you would apply clr to the original compositional data. However it does not follow that clr is applicable to distance data derived from compositional data. (distances are not compositions.) Why do you need to transform the distances at all once you have derived them?

I need to transform the distances to make the PCoA plot. The way the qiime2 workflow is set up you create the distance matrix file from a feature-table.biom file that is zipped as a .qza file. This makes the actual feature table hard to manipulate, and I have been having trouble exporting that, transforming it, then re importing it into qiime2.

https://github.com/LangilleLab/microbiome_helper/wiki/Amplicon-SOP-v2-(qiime2-2022.11)
Here is the workflow I am using in qiime2. The problem is occurring at step 8. As you can see, the deblur_table_final.qza is the input for the diversity metric step (this calculates a whole bunch of diversity tables and graphs, including the beta diversity PcoA plot I want). But as I mentioned before, this file is hard to manipulate & apply the clr() function to. An output of this step is a distance matrix, and I was told on another forum it was possible to transform the distance matrix in R.

I am still not clear on why you need to transform the distances (there is no mention of that in the linked doc?), but if you did you could use some of the suggestions from the qiime2 team here (note these do not include CLR). I would suggest figuring out how to apply transformations to the feature table (for which you should feel free to post questions here of course!), not transforming the distance matrix, or transforming with one of the recommended methods from that post.

Based on conversation in the comments.
CLR is intended for compositional data. Distance data are not compositional, and compositions are (generally) not symmetric. So you should perform your normalization differently.

If you want to transform your data for normalization purposes on the features themselves (in your case, the organismal relative abundances), you should do it on the compositional input from which you generate your distances. This is the intended application for CLR. For a reference on normalization approaches for compositional data (including CLR), see here.
If you want to transform your distance matrix (e.g. the relationships between samples, irrespective of the underlying relative abundances), then you should use one of the recommended methods for transforming distance data.

Note that the applications are very different. Application (1) is transforming the underlying organismal frequency data, such that all of the relative abundances are on the same scale. This will make all downstream analyses more stable / interpretable, including ordination analysis.
Application (2) is saying basically, "my ordination / PCoA on my set of samples still looks weird, even after all that normalization I did on my relative abundances. I want to fix it." Examples of "looks weird" given include negative eigenvalues.
You can't substitute (2) for (1). They are working with different data and have different consequences, though both can lead similar issues in your outputs and ordination plots
It seems like your fundamental problem is that qiime2 outputs are difficult to work with. If that is the case, I would suggest asking a separate question about how to handle qiime2 outputs (or contacting the qiime2 authors about how they do compositional normalization).
It looks like the qiime2 forum already has questions about the right way to apply CLR or related transformations in their workflow. Here is another one. I would suggest following those recommendations and spending some quality time with the tutorials and user forum.


text:
Removing missing categories from geom_bar
I have a simple dataframe that includes points where the x values are 0, 1, 2, 3, 4, and 7. When I go to plot the data this is what ggplot2 gives me:

I only want to include the points mentioned above. How do I take out the instance of 5 and 6 in the bar plot?
(I have tried googling the answer and nothing is coming up with my actual need so I figured I'd ask the community which might give an answer sooner - please be kind I know this is bound to have a simple solution)

I’m voting to close this question because it’s not a bioinformatics question and should be on stack overflow instead.

I was able to get it to work by converting the column in the dataframe to a character using this command:
df$day <- as.character(df$day)
and it worked like a charm! Thanks for your advice @haci

You can give a try by converting the column of numbers to a factor variable.


text:
How do I correctly format my limma/eBayes code?
Here is part of my data table, which I will call "df":
    Subject.ID    Sex           V1           V2           V3           V4           V5
1   GTEX-1117F female   4.24944534   0.18307358   1.99276843   0.21785110   0.74777407
2   GTEX-1128S male    3.286585483  0.165944088  1.843983844  0.444455046  0.307311939
3   GTEX-11EMC female   3.31947343   0.27846061   1.46519089   0.33753984   0.91708799
4   GTEX-11GSP female  3.232353768  0.215492230  1.778208576  0.166458372  0.644225472
5   GTEX-11TTK male    2.934705200  0.256406854  1.712815854  0.256165277  0.948301812
6   GTEX-11ZUS male     3.77188558   0.18434378   2.43189035   0.37373172   1.14339342
7   GTEX-12WSD female   4.22032995   0.21933892   0.93900085   0.08687886   1.06901468
8   GTEX-12ZZX female   3.43616184   0.39473363   2.73205207   0.33525457   1.51399598
9   GTEX-1313W male     3.66334462   0.39418485   2.47248777   0.53545580   0.64634702
10  GTEX-131XW female  2.956614288  0.276317993  1.977096712  0.167743280  1.528071165

How do I correctly format the following code to account for the kind of dataframe I'm working with? I'm using sex as the factors to be interacted. Here is what I have so far:
design <- model.matrix(~ Sex)
fit <- lmFit(df, design)
fit <- eBayes(fit)
topTable(fit)

The first line is giving me the error "Error in eval(predvars, data, env) : object 'Sex' not found" (and idk why I'm getting this error since there is a "Sex" column in df). The second line gives me the error "Expression object should be numeric, instead it is a data.frame with n non-numeric columns" (and I know why I'm getting this error; I'm not sure how to handle the error, though).

You need to look at how to create a model matrix : https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/model.matrix
In your code, you did not specify the data from which the Sex columns should be taken, or you did not created the Sex variable
Either create Sex=df$Sex or specify it in the model.matrix function model.matrix(~ Sex, data=df)

Thanks for the response

I'm guessing you didn't address the other part of the question because this post had 2 questions embedded in it. If you know how to address the other part, could you help with that as well, please? I made another post while taking out the part that you answered.

It is impossible to apply this regression on non-numeric column (error is pretty clear), you need to select the adequate columns

How would I go about doing that? I've tried lmFit(data.matrix(new_dataset[,-c(1,2)]), design), which gives the error ""row dimension of design doesn't match column dimension of data object"" and fit <- lmFit(new_dataset[,-c(1,2)], design[,-c(1,2)]) which gives the same non-numeric column error.

Matrix to regress should be "A matrix-like data object containing log-ratios or log-expression values for a series of arrays, with rows corresponding to genes and columns to samples." You cannot regress the df object as it stands. It should be a matrix with samples as columns and V1, V2, etc as rows. You need to format df accordingly, everything is written in the lmFit description https://rdrr.io/bioc/limma/man/lmFit.html

Okay, I transposed my df (I think that's what you meant). I'm still getting the same numerical error (what I have currently is fit <- lmFit(transposeddf[-(1:2)], design)) (skipping the first two rows which are non-numeric). And the reason the matrix is flipped is because I needed to transpose it earlier for a separate process.


text:
Looking for datasets of RNA-Seq (or scRNA-Seq) for L1 larvae C. elegans
I was wondering if anyone know of a repository for C. elegans data sets, specifically in L1 larvae stage.
Are there any papers or GEO-sets for this organism at this stage?
thanks
Assa

See WormBase, which has a number of different resources for different nematodes:
https://wormbase.org


text:
How to parallelize BEAST
I'm running BEAST v1.10.4 with BEAGLE v4.0.0. My dataset has one large partition and I have access to a machine with many CPUs so I would like to optimize performance by parallelizing the process. Following the documentation and various (quite old and probably outdated) help Q&A pages I tried
beast -beagle -beagle_CPU -beagle_instances X -beagle_thread_count Y -threads Z mydata.xml

with different values for XYZ but BEAST is ever only using a single CPU regardless.
How do I make BEAST use on more than one CPU?
We have a 'server' in-house, that is essentially a big computer with 144 CPUs that is shared by all users. That's why it would be important to limit the max amount of threads with the -threads option. No queuing system, no GPUs. I want to do some tip date sampling which afaik is not available in Beast2.

The machine is important because of the introduction of M1 type chips which will run GPU.
From historical use of beast1 the following will work ...
beast -overwrite -beagle_auto -seed 2001 mydata.xml > myfile.runtime.log

This will parallelise happily within the number of processors available. I would recommend it. This is because it may not parallelise across 144 CPUs. It depends on the machine architecture. Also from memory (mine not RAM) it dynamically readjusts, so even if it was 144 CPUs parallelised when another job lands on a CPU that drops to 143 CPUs for beast1.
Anyway, I use the Cloud so I no longer keep up with computer architecture.
If this floods 144 processors and there's no dynamic readjustment I would suggest,
 beast -beagle -beagle_SSE -beagle_CPU -threads 10 -overwrite mydata.xml > myfile.runtime.log

I would limit the number of threads to get it working before considering e.g. 40.
However, shifting to GPU is strongly recommended, this will solve any beast2 issues and its very fast. However, I don't know about beast1 for GPU, I've never used it for that.

If none of the above work and you've not get any sort of M1 or higher chip (to move to GPU) then you need to look at the error log, especially if you see this error.

Error 122 parsing the xml input file
Cannot create class: ThreadedTreeLikelihood.

Note although you've got a nice CPU heavy machine, simple fairly modern Mac will do the same operation probably quicker via GPU, at least for beast2, beast1 I dunno

Just to note I'm acting in solidarity with an SE strike until an internal matter is resolved about ChatGPT and requested a meeting with staff to explain why. I've formally notified SE. I'll be answering very select questions and doing only essential admin  (breaches of CoC and spam). Hopefully all is back to normal shortly.

Thanks! Unfortunately I couldn't get either to work (the first suggestion doesn't run at all), however, I followed your general advice to use more/different resources to speed up beast worked for me.


text:
Extract data from genome
This question was also asked on Biostars 3 months ago, i.e. a long time - so its perfectly reasonable for the question to be asked here
I started studying bioinformatics and I enjoy it a lot (I'm a software developer with around 15 years of experience but mostly in backend, scripting, architecture, and security) and I want to do, as a project, a database regarding natural compounds. When I saw that on NCBI there are a lot of plant genomes I thought, that, with software, somehow we can "guess" (probabilities, not certainties) what a plant can contain or not.
I'm wondering, if, by bioinformatic means only, it is possible to "extract" genes and, after that proteins and natural compounds from plant genome data that is available on NCBI (please see https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_029618835.1/ as an example)? Example: just by processing the data from the link above we can get to Flavoxanthin or any other substance.
I was wondering if a kind of "de-novo" annotation can be made (if that's the correct term).
What other useful information from a chemical compound perspective can be found from digging on the genome itself?
Thank you!

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

Why do you want to do this? There are lots of different tools; it would help to be more specific in order to make the answers actually useful.

@gringer .. you are totally right, I put full details on another platform when I saw that the interest here was almost none:  https://www.biostars.org/p/9566425/

Thanks for the additional information, but your additional details still describe a generic problem, and are therefore unlikely to get good answers. What is the specific problem that you are trying to solve?

This sounds like you are wanting to perform a machine learning calculation. However, for this question leveraging KEGG, Kyoto Encyclopaedia of Genome and Gene and building an SQL database around it is much more effective. It represents around 60% of all functional genes. KEGG is the best gene annotation out there IMO, particularly for flavonoids etc... Fantastic and all their linking pathways, everything links to everything, just amazing. KEGG represents an enormous volume of work, much of it being manually crunching publications for each gene.
The entire database of KEGG keys is downloadable. Connecting a gene to a class of function is easy via KEGG the relationships are defined its just SQL-ing the whole thing together (querying all the tables). Again, it's totally and completely brilliant. Sure they'll be stuff out there has SQL-ed KEGG key tables, but if your a developer its better to build it independently - this way you're certain whats in there and what its doing.
However, linking raw sequence data to a KEGG gene code is more challenging. lastKOALA and GhostKOALA is ... yeah well an approach no more. There's a main stream approach but with current limitations in implementation. KEGG however isn't limiting, it's linking a gene to a KEGG id is the challenge. Getting hold of the raw sequence for KEGG would be great, fantastic in fact, but I'm yet to find it.

I've used Augustus stand alone (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1160219/) to find genes in a given eukaryotic genome. I agree you can then use things like KofamSCAN to categorize the genes against the KEGG database. If you want to compare against the NCBI database, diamond blast (https://github.com/bbuchfink/diamond) is probably the quickest way to go :)


text:
Nextflow: pass reads variable to process in another script
This question has been asked and answered on Stack Overflow
I'm new to nextflow, channel, process. Kindly bear with me. I'd like to create a pipeline, pass paired reads from one to another process. I am unable to do it.
script1.nf
params.outdir = 'results'

params.reads1 = "/Users/user/Downloads/tiny/normal/*_R1_xxx.fastq.gz"

params.reads2 = "/Users/user/Downloads/tiny/normal/*_R2_xxx.fastq.gz"

println "reads: $params.reads1"

println "reads: $params.reads2"

Channel.fromPath(params.reads1,checkIfExists: true).view()

Channel.fromPath(params.reads2,checkIfExists: true).view()

include { fastp } from '/Users/name/Documents/name/nextflow_scripts/fastp.nf'

fqpairs_ch=channel.fromFilePairs('/Users/name/Downloads/tiny/normal/*_R{1,2}_xxx.fastq.gz')

process check {

    output:
    stdout

    """
    echo "hello from command line \n"

    """

}

workflow {

    check().view(). /// this works

    fastp(fqpairs_ch).view() // it gives error

}

fastp.nf
params.outdir = "./results"

process fastp() {

tag "fastp on $sample_id"

input:
        tuple val(sample_id), file(x), file(x1) from fqpairs_ch

output:
println ("$sample_id")

    script:

    """

mkdir fastp_trimmed

        fastp \
        --correction \
    --cut_tail \
    --disable_trim_poly_g \
    --length_required 50 \
    --qualified_quality_phred 20 \
    --thread 12 \
    --trim_poly_x \
    --unqualified_percent_limit 20 \
        -i ${x[0]} -I ${x[1]} \
        -o fastp_trimmed/joinedfiles.dat_${x[0]} -O fastp_trimmed/joinedfiles.dat_${x[1]}

    """

}

Error:

ERROR ~ No such variable: fqpairs_ch
-- Check script 'fastp.nf' at line: 13 or see '.nextflow.log' file
for more details

I'm unable to wrap by head around tag, val. Also, it is difficult for me to understand how to pass variable from one script to another. I've watched Youtube tutorial but with limited success.


text:
The mathematics of kinship (degrees and distributions of relatedness on individual and group levels)
I have been thinking about to what extent a biological entity’s characteristics can be decomposed into discrete sections of the genome that are parsed and mapped one-to-one to a physical structure. It is fascinating in this regard how indistinguishable biology is from computational linguistics and formal systems theory in general - DNA is a formal language; and we are still trying to unravel the grammar. Right?
I think the following question depends on what type of “grammar” DNA has. I don’t know the term in the field, but I have been using the term “holographic” in my own thinking, when a sequence in a code is interpreted twice - by two different grammars, or two different parts or stages of a grammar, contributing meaning to multiple realizations in the output. (“Holographic”, as in, a 2D surface appears to be 3D so long as the information is encoded on it in such a way to be parsible / reconstructible - the information is not lost, in the mix, but is able to be separated and mathematically resolved (some sort of non-associativity under some function, is what I suspect).
My question is about to what extent “kinship” is indicated by “having similar genes”.
Genes are virtual (I think) - pure information. It seems like the same genes can come into being without them having replicated from the same ancestor (convergent evolution). The presence of the same gene in two biological entities does not strictly imply closer “kinship” when compared to two biological entities without that shared gene. Conversely, we need not expect entities of close kinship to share the same genes (a sibling can inherit one allele of some gene, another sibling the other one, dominant or recessive, etc.)
Perhaps there is an analogy between language and biology, that I haven’t clearly spelled out, yet.
-> Organisms are to DNA what sentences are to “logical form” (in linguistics).
-> DNA is to mechanisms that build DNA (RNA transcription, protein synthesis, etc.) what logical form is to “universal grammar” (or, the actual system cradling the human aptitude for languages, itself).
-> (The odd part): Biological reproduction is when a “sentence” of language directly, physically transforms (or, contributes a part of itself) into an instance of a “grammar” (a valid grammar, as specified by the parameters/constraints/dictates of “universal grammar”). It is intriguing how the analogy gets challenged here, perhaps to be probed further. For now, it could be seen to stretch to how human language can express grammars (like humans can define grammars - constructed languages, as one example - mathematics - or any formal system at all.)
-> Lastly, biological evolution has some analogy to language evolution (perhaps, a fourth “layer” or “grammar” (a specification of means by which something can change - rules, constraints) of the system), in that there is some stochastic generation which gets selected for some (chaotically changing, varying) context (ie, natural selection; as for humans and language, it might be hard to specify the various reasons, including highly social ones, why language features, or cultural features, are either retained or dissipate (Dawkins’s concept of “memes”)).
I think there are two competing dynamics.
One says that similar genes are not “informationally significant” because the “universal grammar” and the transformation rules/rounds/generations on universal grammars themselves (the ways the system can vary from generation to generation, in the space of possible evolutionary trajectories) can easily generate the same structure with multiplicity, independently.
Another competing force says that presence of similar genes indicates a likelihood of kinship because it is improbable that two people who collectively have a high degree of similarity between their total collection of genes, perhaps (this is how ethnicities arise - a group exchanging genes internally, within a single evolutionary path, pass the same genes around amongst each other; the same with human languages).
I wonder if there is any “equilibrium” theory, where the extent to which a system’s propensity to permit random duplicates of a form, is somehow balanced with the (evolutionary) dynamic where memes (genes, language features) cluster with regards to proximity of interaction.
So we at the end of the day have a situation where similar overall gene profiles is actually a probable indicator of “poly”-kinship (a kinship vector where any two entities have some last common ancestor for each of their genes, not just their last common “parent” - and those kinship vectors can have unique distributions, from homogeneous to with wide range - but the same average).
But, there are still outliers where a person is able to have a similar phenotype because of sufficient randomness in the way genes work, without it indicating a high-value (normed) kinship-vector.
(This could apply to linguistics where it would be fallacious to assume Japanese and Turkish must be related to each other just because they share strikingly agglutinative morphology features.)
How is this modeled mathematically, in bioinformatics?
Thank you.

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

The question is too broad and needs to be focussed. Game theory has been applied to population dynamics to explain diverse behaviours of "fitness". Phylogenetics is the tool to answer to your question and has been applied to linguistics.

By game theory I mean understanding equilibrium selection.

Could you please suggest good sub-questions the post could be broken down into? Experienced people in a field often are able to ask better questions because they have more meta-knowledge, but a beginner has many strands of interconnected thought they can’t organize enough yet to separate. I guess my question is about the “formal grammar” of biology (to what extent is it known). It is quite broad. I guess I’ll polish off new sub-questions, one at a time. Yes, that could be a good approach.

Do you have a specific project or analysis goal in mind? Because bioinformatics is such a broad area of study, It's usually better to start off with your actual problem, and find tools or solutions to fix that.


text:
Extract compounds/proteins from plant genomic data
I want a process by which I can identify genes within plant genomes; i.e. a process that works for all plants, not just for a single species.
I'm wondering, if, by bioinformatic means only, it is possible to "extract" genes and, after that proteins and natural compounds from plant genome data that is available on NCBI (please see https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_029618835.1/ as an example)?
Example: just by processing the data from the link above we can get to Flavoxanthin or any other substance.
What other useful information from a chemical compound perspective can be found from digging into the genome itself?
Any hints will be greatly appreciated.
Thank you!

I believe I understand the project, to obtain a biochemical/metabolic profile of a large number of poorly annotated plant genomes. There isn't a perfect solution, like Prokka for bacteria. It is clearly an area subject to recent development and interest.
I would recommend using Gene Ontology Meta Annotator for Plants (GOMAP) which is described in "Plant Methods" journal.
There are limitations:

KEGG is not used (below);
its a singularity container - and thats its only distribution.

Learning Singularlity You will need to learn running the container through Singularity and that will be limiting if you are not familiar with containers. However, this is the best of the available "out of the box" solutions I can see in plants. A Singularity HPC is a good call here, tricky if you don't have access to one. I know Docker, Singularity is similar but different.  The information on Singularity is here. I would start reading from "Executing Commands"
singularity pull docker://godlovedc/lolcow # get the container
singularity exec lolcow_latest.sif cowsay moo # run it

(the cow is a famous teaching tool for containers)
KEGG One approach that can be used is to perform COG (NCBI/Genbank code) and KEGG no. identification (e.g. eggNOG). Tools such as eggNOG are setup for metagenomics data however. The COG/KEGG no. output is only the starting point, interrogating each of these to assign gene function and metabolic pathways will require a KEGG database and custom code. Unless someone has published, e.g. KEGG database in SQL.
In my experience the Kyoto Encyclopaedia of Genes and Genomes (KEGG) is the best singular functional database available, thats why I see GOMAPs exclusion of this as a limitation. COG doesn't provide a particularly useful information functional description IMO. eggNOG, or any other method to identify COGs  and KEGG numbers, interrogated against the KEGG database works for plants very well. I know this because I'm not focused on plants but often recover plant  defence mechanisms genes from KEGG. Basically, KEGG identifies >> half of all genes with known COG nomenclature and notably higher than other functional databases.Thus I see GOMAPs exclusion of KEGG as a limitation. However, to reiterate to go identify COG/KEGG -> KEGG database interrogation and analysis is a lot of custom code.
Databases GOMAP implements Arabidopsis annotation, which may be cool and provides a starting point. It's also running UniProt, UniProt-plants, interpro, GO databases and those + Arabidopsis annotated DB could be really cool. It definitely sounds promising.
Processor time Any method of gene identification takes AGES and I have no doubt that GOMAP will be REALLY slow, looking at the pipeline it's using - there are a lot of databases its using. This is GOMAPs justification for using Singularity, it going to take ages to run therefore containerise it in a format often used by modern HPCs. To be honest it's a bit lazy, especially as it was published in 2021 (should be uploaded e.g. to conda). The authors talk about Python and R ... hmmmm, it looks to be a MySQL backend, thats where its speed will come from. The container will be large BTW hopefully is <1GB, but it will be >200MB.
Summary Its worth learning Singularity just keep to the basic commands to run GOMAP. You will need access to a large machine and the output will be okay but they're better general databases. However, I don't know how useful the Arabidopsis database it - thats about the biology, i.e. what you are looking for and whether Arabidopsis is a good model of that investigation
Note, Springer Nature refused a link image request for the container overview, but the image I wanted to share is here.

Just to mention the essential missing part the Singularity container test described singularity exec lolcow_latest.sif cowsay moo. The output is moo. Cows say "moo".

Thank you for pointing to GOMAP, seems a good fit for my plan. Basically, I'm doing a plant database (already have an online demo) as my Master's final project and I wanted to add these detectors besides aggregating other data.
I never heard about Singularity before, but, at base, I'm a software developer and I worked with docker, vagrant, virtualbox and other similar systems so it shouldn't be an issue to get it rolling.

Thats good. The project makes a lot of sense. Singularity is huge in biological science (Docker uptake was low). The only thing to add is entering the Singular container is very different to Docker's -it. From there both the code and databases are open for inspection.


text:
What statistical analysis should I perform on DNA variants from unequal number of case (34) and control (11) samples?
I have performed NGS on 34 patient samples and 11 controls. Now I have variants from both groups but there is obvious difference in number of variants.
How can I compare variants of both groups statistically?
Can I do association analysis?

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 simply a GWAS calculation. You can do this via phylogenetics, it's very complicated, but it will work. HOWEVER, I am assuming this is a whole genome NGS (variants so it must be), rather than RNAseq. RNAseq is a very different calculation. It would be useful to specify the disease because there will be prior variant association studies.

Please explain more about your experiment and desired analysis. What platform was sequencing done on? What coverage? How do you know you have variants in both groups? What additional processing has been done on the FASTQ files to demonstrate that there are variants?

Mitochondrial DNA sequencing on Illumina platform and analysis was performed using bwa-mem, bcftools and Mitomap for variant annotation. I want to compare DNA variants in both groups (34 patients, 11 controls) and any association of these variants with disease. Can population specific Allele frequency from gnomAD help in analysis of rare variants?


text:
How to identify index sequences for cutadapt for atac-seq
I am just starting on my atac-seq analysis (not very experienced, so my apologies for a stupid question), and I am wanting to trim the adapters of off my paired-end reads.
I am following the galaxy atac-seq protocol that says to enter in the custom 3' adapter sequence. I took the primer sequences from Buenrostro.
>Ad1_noMX
AATGATACGGCGACCACCGAGATCTACACTCGTCGGCAGCGTCAGATGTG
>Ad2.7_CTCTCTAC
CAAGCAGAAGACGGCATACGAGATGTAGAGAGGTCTCGTGGGCTCGGAGATGT

Following the ATAC protocol, I see that the primer sequences they used are the reverse complement sequences of the end of the provided sequences (ACATCTCCGAGCCCACGAGAC/CACATCTGACGCTGCCGACGA). I took this part and added the reverse complement index primer of 2.7 (GTAGAGAG) to the front.
So for my R1 my sequence is:
GTAGAGAGACATCTCCGAGCCCACGAGAC
And for my R2 my sequence is:
GTAGAGAGCACATCTGACGCTGCCGACGA
However, when I run the command I only get 1.4% read1/2 with adapter. Could anyone suggest where I am going wrong?
Thank you

The sequence to trim for ATAC-seq is CTGTCTCTTATACACATCT. It's the Nextera adapter sequence. The sequence is the reverse complement of the 3' ends of your primers, both forward and reverse.
You trim adapters which come before any index sequence so trimming adapter (and by this everything beyond it) is sufficient.


text:
How are the "effect_allele" and "other_allele" bases chosen in a PGS file from the PGS Catalog?
How are the "effect_allele" and "other_allele" bases chosen in a PGS file from the PGS Catalog?
For example, for PGS002723 (https://www.pgscatalog.org/score/PGS002723/), if I set the "effect_allele" to "ref", then "ref" is equal to the reference
base of hg38 in 98.60% of the cases on chromosome 22.
On chromosome 22, there are 205 PGS variants where the "effect_allele" does not match the
reference base of hg38. In all of these 205 cases, the "other_allele" matches the hg38
base.
Why does PGS Catalog encode the hg38 reference base in 14,477 of cases for chromosome 22
as "effect_allele" and in 205 of cases as "other_allele"?
Finally, how can I deal with this encoding if I want to compute the polygenic score on a certain cohort?
Can I simply invert the sign of the "effect_weight" column for these 205 cases where the "other_allele" matches the reference?
I paste my reproducible code below:
suppressMessages(library(BSgenome.Hsapiens.UCSC.hg38))
suppressMessages(library(data.table))
suppressMessages(library(R.utils))
suppressMessages(library(stringr))

setwd('/cluster/home/criccio/PRS/')

source('scripts/get_hg38_base.R')

if (!dir.exists('output/table'))  dir.create('output/table')

PGS_ID <- 'PGS002723'

# Example to ask questions online
system(paste0('wget -P output/ --recursive --no-parent ftp://ftp.ebi.ac.uk/pub/databases/spot/pgs/scores/', PGS_ID, '/ScoringFiles/'))

scoring_file <- data.table::fread(paste0('output/ftp.ebi.ac.uk/pub/databases/spot/pgs/scores/', PGS_ID, '/ScoringFiles/Harmonized/', PGS_ID, '_hmPOS_GRCh38.txt.gz'))
data.table::setnames(scoring_file, c('chr_name', 'hm_pos', 'effect_allele', 'other_allele'), c('chr', 'pos', 'ref', 'alt'))

pgs_variants_per_chrom <- scoring_file[, .N, chr]
data.table::setnames(pgs_variants_per_chrom, 'N', 'nb_pgs_in_scoring_file')

chromosome_number <- 22

scoring_file[chr == chromosome_number, ref_hg38 := sapply(scoring_file[chr == chromosome_number, pos], get_hg38_chr22_base)]
scoring_file[chr == chromosome_number, ref_hg38 := toupper(ref_hg38)]
scoring_file[chr == chromosome_number, table(ref == ref_hg38)] # 14477 TRUE and 205 FALSE
scoring_file[chr == chromosome_number, 100 * mean(ref == ref_hg38, na.rm = TRUE)] # 98.60%

scoring_file_mismatch <- scoring_file[chr == chromosome_number & ref != ref_hg38]
scoring_file_mismatch[chr == chromosome_number, mean(alt == ref_hg38)] # 100% match
data.table::fwrite(scoring_file_mismatch, paste0('output/table/', PGS_ID, '_chr_pos_ref_mismatching_hg38_chr', chromosome_number, '.csv'), sep = ',')

get_hg38_base.R
library(BSgenome.Hsapiens.UCSC.hg38)
Hsapiens_chr22 <- strsplit(as.character(Hsapiens$chr22), '')[[1]]
get_hg38_chr22_base <- function(pos) Hsapiens_chr22[pos]

Effect allele is almost always the reference allele. It theoretically does not have to be, so it's good to not rely on that fact. However, in practice it is typically a safe bet. At first glance, the answer to your question could be this:
The polygenic score from the PGS catalog is using GRCh37. Your genome build is hg38. I expect that some alleles in the score swapped ref/alt positions between builds, resulting in your 205 cases where the effect allele matches the alt allele.
To confirm, could you list a few examples of cases that matched as expected and a few examples of cases that didn't match?
You don't need to deal with this encoding if you want to compute a polygenic score on a cohort. You can match based on the allele itself. (Also, there are existing tools which will calculate polygenic scores for you.) Here is an example (just a snippet, you'd have to iterate through your files, etc. for this to work):
            if PGS_row scoring file identifier like rsID or position == vcf_row identifier: # When the SNP is in both files
            effect_allele = PGS_row[2] # Wherever the affect allele is in your scoring file
            effect_weight = float(PGS_row[4])
            dosage (of ref allele) = (depends on your VCF formatting)

            if effect_allele == vcf_row[3]: # Effect allele matches reference allele
                running_total += (2 - dosage) * effect_weight # We assume dosage of alt is 2 minus dosage of ref. So if you have 2 alt alleles, 2-2 = 0 and there will be nothing added for this SNP.
                num_additions += 1
            elif effect_allele == vcf_row[4]: # Effect allele matches alternate allele
                running_total += dosage * effect_weight
                num_additions += 1
            else:
                nomatch = nomatch + 1

# Calculate final output score
final_output_score = running_total / num_additions
percent_matched = (num_additions)/((nomatch)+(num_additions))

You shouldn't invert the sign of the "effect_weight" column for these 205 cases where the "other_allele" matches the reference. Not having the effect allele doesn't add the opposite beta to the score, it simply doesn't change it.
Try this:

Do you still see those 205 cases if you use a genome build matching the PGS file? (I.e., GRCh37.)
Do you still see those 205 cases if you use the LiftOver PGS scoring file with hg38 coordinates? (Download it here.)

Thanks a lot. This answer is really helpful. Yes, I checked the scoring file of GRCh37 and there were no mismatches, so it must be as you say that the reference base changed between GRCh37 and GRCh38. I also thought inverting the sign was not OK but now you convinced me. Indeed, it add weight if you have the effect allele. Cheers!


text:
datasets or biobanks with both (WGS or genotype intensity) and methylation data
My aim is to associate CNVs to DNAm, but I don't know where to find datasets which contain both (sequencing + methylation). Do you know any?


text:
consortium, biobank or browser with CNV eQTLs
I found that GTEx contains several SNPs eQTLs, but I'm looking for previously reported CNV-eQTLs, where can I find such data (freely available)?


text:
How do I fix chromosome-specific Michigan Imputation Server error (Imputation on chromosome 14 failed. Imputation was stopped)
I am imputing human genome samples, all hg19, using the Michigan Imputation Server. The input of the data is one VCF per autosome, which each VCF containing a few dozen samples.
The samples within the VCF are from highly heterogenous sources (various genotype arrays, from multiple labs). I have cleaned and standardized the data (removed indels, multiallelic sites, duplicates, and non-standard genotypes). I merged the samples into a single file for each chromosome. The VCF may not be left-aligned.
It passes all QC with no issues, and then fails and gives no error message beyond "unfortunately, your job failed" and "Imputation on chromosome 14 failed. Imputation was stopped."
As a test, I imputed chromosome 1, chromosome 13, and chromosome 14  by themselves. Chromosome 1 and 13 worked perfectly, but chromosome 14 failed imputation. I have no idea why chromosome 14 would be special or different. Also, imputing chromosomes 1-13 failed on chromosome 9, so it appears only certain chromosomes are causing issues.
Does anyone have ideas as to why this might be occurring, or how to fix this issue?
I dislike using hg19. I am using that because that is the build that some of the genome data I received is on. Labs are still using hg19 in the present day, for some reason. However, I have used hg19 successfully on the Michigan Imputation Server before, and it supports both hg19 and hg38, so I don't think build is the issue. I have a different batch of hg38 data that I will run through the server later.
Here is a comparison between the QC statistics:

Chromosome:
13 (success)
14 (failure)

SNPs
75,675
69,049

Chunks
6
6

Datatype
unphased
unphased

Build
hg19
hg19

Reference Panel
apps@1000g-phase-3-v5 (hg19)
apps@1000g-phase-3-v5 (hg19)

Alternative allele frequency
> 0.5 sites: 14,015
> 0.5 sites: 11,776

Reference Overlap
99.18%
99.25%

Match
44,261
39,947

Allele switch
0
0

Strand flip
0
0

Strand flip and allele switch
0
0

A/T, C/G genotypes
0
0

Filter flag set
367
342

Invalid alleles
20,567
18,022

Multiallelic sites
12
16

Duplicated sites
0
0

NonSNP sites
0
0

Monomorphic sites
9,963
10,302

Allele mismatch
138
116

SNPs call rate < 90%
43,096
38,760

Excluded sites in total
74,143
67,558

Remaining sites in total
1,165
1,187

Typed only sites
367
304

Execution status
successful
successful

Per the Michigan Imputation Server documentation, I've checked my input VCF files with CheckVCF.
Here are the results on the chromosome VCF file that failed:
---------------     REPORT     ---------------
Total [ 69095 ] lines processed
Examine [ 49 ] VCF header lines, [ 69046 ] variant sites, [ 21 ] samples
[ 0 ] duplicated sites
[ 0 ] NonSNP site are outputted to [ test.check.nonSnp ]
[ 0 ] Inconsistent reference sites are outputted to [ test.check.ref ]
[ 0 ] Variant sites with invalid genotypes are outputted to [ test.check.geno ]
[ 11947 ] Alternative allele frequency > 0.5 sites are outputted to [ test.check.af ]
[ 34408 ] Monomorphic sites are outputted to [ test.check.mono ]
---------------     ACTION ITEM     ---------------
* No error found by checkVCF.py, thank you for cleanning VCF file.

Here are the results on the chromosome VCF file that succeeded:
---------------     REPORT     ---------------
Total [ 75718 ] lines processed
Examine [ 49 ] VCF header lines, [ 75669 ] variant sites, [ 21 ] samples
[ 0 ] duplicated sites
[ 1 ] NonSNP site are outputted to [ test.check.nonSnp ]
[ 0 ] Inconsistent reference sites are outputted to [ test.check.ref ]
[ 0 ] Variant sites with invalid genotypes are outputted to [ test.check.geno ]
[ 14059 ] Alternative allele frequency > 0.5 sites are outputted to [ test.check.af ]
[ 37553 ] Monomorphic sites are outputted to [ test.check.mono ]
---------------     ACTION ITEM     ---------------
* No error found by checkVCF.py, thank you for cleanning VCF file.

I would suggest chunking chromosome 14 and performing each chunk sequentially. I don't do eukaryotes full stop, but thats how I'd do it. I think there's a getter / setter absent in an algorithm.

Interesting. Why do you think chunking might work? And I'm not sure how to handle that with imputation. Perhaps it impute "beyond" each chunk somehow, and then combining again would be difficult?

Generic divide and conquer (how coders think), breaking the problem into smaller pieces ('chunks'). If its a data issue, one 'chunk' would fail. If it was RAM (unlikely) chunking would get around the problem. ... Again I don't do eukaryotes, esp. the human genome. Someone must have done this already via bcftools, vcftools etc ... the bits of the human genome I've done I noted I could code it just as easily as existing tools.

Okay, got it. Good idea; chunking would help narrow down where the problem is in the data.

Yep thats it exactly and then compare 'chunks' which work, which don't.

If you're chunking, make sure to create overlapping chunks. The resulting haplotypes can then be merged together by their overlapping regions using standard tools. I wrote one a while back (https://bitbucket.org/wkretzsch/hapfuse/src/master/), but bcftools should do the trick as well.

@BigMistake I've deleted your cross-post on Biostars because you have this bad habit of asking every question on both places.

@RamRS is that against the rules? I looked and couldn't find guidelines, and also waited until the next day, to post there as there are no answers here. If you have a page that details that / other rules, let me know so I don't break any of them in the future.

It's not in the rules, it's common etiquette. Also, you were asked to not do that once, so if you had a problem with that, you should have addressed it then. There is a meta post on this BTW: https://bioinformatics.meta.stackexchange.com/questions/265/about-posting-queries-in-stack-bioinformatics

@BigMistake Stackexchange will not act against cross-posts because that requires a change in the  Code of Conduct (CoC) which is determined centrally. Biostars belongs only to Biostars and has autonomy to do what it wants. However, here we will put a cross-post link into the post. Thats fair for all concerned. I agree cross-posting isn't ettiequte but it is common.

@M__ Thanks. I was having a bit of trouble understanding why cross-posting wasn't liked. I think I mostly get it now. I will make sure to add urls. Though personally, adding content to both sites seems like it would make it more searchable and be an improvement, as long as the answer is copied to both sites with credit to the original author of the answer, so maybe I don't quite understand the ettiequte.

@BigMistake I'm afraid you can't copy and paste answers between sites. The  answer belongs to the author via SE's copyright, so in theory only the author could do that - but no-one else. I don't know how Biostars copyright works however. E.g. if the author posts here first and then copy-pastes to Biostars - thats fine. However, if Biostars own the copyright not the author that author could  not copy-paste their answer from Biostars to here.

Good point. Then, saying "answer here" with a link to the original site at least in a comment might be better.

@winni2k When I try chunking, 100% of the chunks work despite the overall file not working. However, I'm not sure if creating overlapping chunks is necessary since I'm using the Michigan Imputation Server--it looks like they combine everything into chromosomes no matter what.


text:
Nextflow - Process has already been used
I am trying to write a nextflow script below:
/*
 * pipeline input parameters
 */
params.reads = "/path/to/fq/files/"
params.guide_library = "/path/to/guidelibrary/"
params.outdir = "results"

log.info """\
    DUAL GUIDE CRISPR QC - N F   P I P E L I N E
    ===================================
    guide library: ${params.guide_library}
reads        : ${params.reads}
    outdir       : ${params.outdir}
    """
    .stripIndent(true)

/*
 * define the `index` process that creates a binary index
 * given the transcriptome file
 */

process READCOUNT {
    input:
    path guide_library
    tuple val(sample_name) path(reads)

    output:
    path "$sample_name"

    script:
    """
    exactcounts.py -inputlibrary $guide_library -read1 ${reads[0]} -read2 ${reads[1]} -prefix $sample_name
    """
}

workflow {
    Channel
        .fromFilePairs(params.reads, checkIfExists: true)
        .set { read_pairs_ch }
   count_ch = READCOUNT(params.guide_library,read_pairs_ch)

}

But am getting the following error:
Process 'READCOUNT' has been already used -- If you need to reuse the same component, include it with a different name or include it in a different workflow context

I believe adding a comma after val(sample_name) on line tuple val(sample_name) path(reads) will solve the issue.
EDIT: Also documented here: https://github.com/nextflow-io/nextflow/discussions/2714


text:
Running Genome Scale Metabolic Model but with different concentration of medium composition
I want to run a simulation with COBRApy for my bacterial Genome Scale Metabolic Model (GEM). I want it to run on different concentration of medium composition such as 0% NaCl, 10% NaCl, and 20% NaCl.
How to implement this type of simulation in COBRApy?
Thanks in advance.

https://cobrapy.readthedocs.io/en/latest/ I don't think you'll find any other shortcuts. Looks easy enough to implement and its cools that its Python.


text:
Cell Ranger returns too few cells!
I recently performed single cell RNA seq using core services and the web summary output from 10x that they have given shows that there are very few cells detected (100-300) cells per sample. Once this same issue happened in 1 out of 3 samples and this time it happened with 4 out of 4 samples.
What might be the issue? The core says the cell preparation was not good but the cells were sorted for viability and then used to make the library. Any help?
Full Summary snap:

I have uploaded the full html file here : https://drive.google.com/file/d/11eQduL4OC4kjC0BfhGp3Nhoci69bRAy6/view?usp=sharing

For one sample, Cell ranger 7.0.1 gave 284 cells while 6.8.2 gave 2200 cells.

Very odd, can you say something about the say something about the sample type?

Also cellranger 6.8 does not exist to my knowledge ;) can you check that

I suggest runnnjng fastqc and maybe some Samtools modules then multiqc to aid diagnosis


text:
How to get full alignment output from BLASTP?
I am running a blastp search:
time ~/bin/ncbi-blast-2.14.0+-src/c++/ReleaseMT/bin/blastp -db nr -query X.fasta -out nr.X -num_descriptions 4000 -remote

but I have 2 problems with the output:
I. Under Sequences producing significant alignments:,
the output is truncated, and I can't read the exact species name:
XP_002906761.1 elongation factor 3, putative [Phytophthora infest...  848     0.0
KAG3108573.1 Elongation factor 3 [Phytophthora idaei] >gb|KAG3151...  845     0.0
KAG2772535.1 Elongation factor 3 [Phytophthora cactorum] >gb|KAG2...  845     0.0
KAG6976665.1 hypothetical protein JG688_00001147 [Phytophthora al...  844     0.0

II.
The above sequences do not show up on the detailed sequence alignments below.
I have tried numerous variations of "outfmt" variations, but now am getting warnings that I'm using the system too much.  So far, output format 0 is the best.
How can I get the full names and not the truncations?
Why aren't the above species showing up later on, and how can I get the full alignment to show?

The problem I had was solved by num_alignments, which is more complete than the first section.  I do not think that the text truncation is solvable.  However, parsing of the latter half solved the problem.
the correct way to get the most complete output is:
time ~/bin/ncbi-blast-2.14.0+-src/c++/ReleaseMT/bin/blastp -db nr -query X.fasta -out nr.X.0 -outfmt 0 -num_descriptions 4000 -remote -num_alignments 4000

Hmmm, firstly can you attempt multi-threading on a remote command via  -num_threads 8...
blastp -db nr -query X.fasta -out nr.X -num_threads 8 -num_descriptions 4000 -remote

It's a worth a try. Otherwise it could likely be a function of the network access or remote server resources. I don't think this is a function of the local RAM. Therefore performing a remote blast to download relevant sections of the NCBI database, convert that into a local database and then perform local blast seems the best strategy.
Multi-threading and -outfmt work well with local blast.
blastp -query query.fa -db test.db -outfmt 6 -max_target_seqs 1000 -num_threads 8 -out blasthits.txt

The alignment is a function of the output formats in the above example its -outfmt 6 that is not useful, but 1 to 4 are relevant as follows.
 -outfmt <String>
   alignment view options:
     0 = pairwise,
     1 = query-anchored showing identities,
     2 = query-anchored no identities,
     3 = flat query-anchored, show identities,
     4 = flat query-anchored, no identities,
     5 = XML Blast output,
     6 = tabular,
     7 = tabular with comment lines,
     8 = Text ASN.1,
     9 = Binary ASN.1,
    10 = Comma-separated values,
    11 = BLAST archive format (ASN.1)

Multi-threading/parallelisation is the way to overcome memory issues, but really effective only in local operations not remote operations. Again, I think there's too much traffic/ or too many resources requested on NCBI server - and thats where the memory error is occurring.

There's no RAM issue, NCBI said they would deprioritize me. I just need to get all of the alignments showing in the 2nd part of the output file, and I'm not sure how to write the command that correctly. So far, option 0 is best, but I've yet to find a way of telling blastp to not truncate the output, and output the alignments for all

I did say I thought local RAM was unlikely to be RAM. NCBI have deprioritised you? Wow.  Biopython does use scheduling so maybe their blast package may be better. -outfmt 1 is probably the best.  I've never encounter truncation in blast so dunno. What I'm recommending is a two stage process and perform the blast locally.


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

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

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

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

    return (10.0 / self.length) * score

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

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

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

MALWMRLLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGFFYTPKTRREAEDLQVGQVELGG GPGAGSLQPLALEGSLQKRGIVEQCCTSICSLYQLENYCN

and print(repr(entry)) is given as

'MALWMRLLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGFFYTPKTRREAEDLQVGQVELGG\nGPGAGSLQPLALEGSLQKRGIVEQCCTSICSLYQLENYCN'

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

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

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


text:
How to join multiple large .tsv (> 20 GB, bed file) files by two identifier?
This question has also been asked on Biostars
I wonder how can I join multiple large .tsv files (over 20 GB each) by using two identifiers (columns) which I used to do using the full_join function in dplyr package for comparatively smaller files. All files are in bed format, where ref_f.bed will be mapped by the other files by col1 and col2.
#merging the files
dflist <- list(ref_f, file_1, file_2, file_3,....file_n)

Test_out <- dflist %>% reduce(full_join, by=c("chr", "start"))

head ref_f.bed
chr1 1 100 .
chr1 101 200 .
chr1 201 300 .
chr1 301 400 .
chr1 401 500 .
chr1 501 600 .
chr1 601 700 .
chr1 701 800 .
chr1 801 900 .

head file_1.bed
chr1 10101 10200 10158 10177 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10158 10177 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10164 10183 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10164 10183 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10175 10194 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10175 10194 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10175 10194 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10181 10200 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10181 10200 Scerevisiae-UniPROBE-Tbf1.UP00299

head file_2.bed
chr1 10101 10200 10158 10177 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10158 10177 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10164 10183 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10164 10183 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10175 10194 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10175 10194 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10175 10194 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10181 10200 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10181 10200 Scerevisiae-UniPROBE-Tbf1.UP00299

head file_3.bed
chr1 29101 29200 29111 29124 Hsapiens-HOCOMOCOv10-AP2D_HUMAN.H10MO.D
chr1 29101 29200 29132 29152 Scerevisiae-UniPROBE-Cha4.UP00300
chr1 29101 29200 29116 29133 Scerevisiae-UniPROBE-Ecm22.UP00305
chr1 29101 29200 29126 29143 Scerevisiae-UniPROBE-Hal9.UP00315
chr1 29101 29200 29126 29143 Scerevisiae-UniPROBE-Hal9.UP00315
chr1 29101 29200 29106 29125 Scerevisiae-UniPROBE-Rpn4.UP00272
chr1 29101 29200 29110 29130 Scerevisiae-UniPROBE-Rsc30.UP00279
chr1 29101 29200 29112 29132 Scerevisiae-UniPROBE-Rsc30.UP00279
chr1 29101 29200 29133 29152 Scerevisiae-UniPROBE-Stp2.UP00338

R needs to import all files before performing the task because the file size for more than 20 files surpasses 20 GB, making this task appear difficult for R. This is the first time I have encountered this kind of issue, and I previously used this method in R for merging smaller size files. Therefore, I'm seeking a better solution to this issue.
It would be a great assistance if anyone kindly advised how to overcome this issue.
Thanks

Can you please give more information about what you're trying to do with these join functions? It's likely that there's a better way to solve your problem that doesn't involve the intermediate step that you're asking about.

If you are open to a GNU/Linux solution, the following should be a good start:
awk '{printf("%s-%s %s %s %s\n", $1, $2, $1, $2, $3, $4);}' ref_f.bed | LC_ALL=C sort -k 1 > sorted_modified_ref_f.txt
awk '{printf("%s-%s %s %s %s %s %s %s\n", $1, $2, $1, $2, $3, $4, $5, $6);}' file_1.bed | LC_ALL=C sort -k 1 > sorted_modified_file_1.txt
join sorted_modified_ref_f.txt sorted_modified_file_1.txt

The awk commands are creating a "first column" that is a combination of columns 1 and 2 and this is required because join "joins" on a single column.

For joining / merging / intersecting bed files, you're better off using bedtools.
It's unclear what you're trying to do with these steps, but I'm guessing it's something like the following:

The bedtools sub-commands include:
[ Genome arithmetic ]
    intersect     Find overlapping intervals in various ways.
    window        Find overlapping intervals within a window around an > interval.
    ...
    merge         Combine overlapping/nearby intervals into a single interval.

If you're interested in combining multiple lines of different files in R, have a look at the map_dfr function (reference here), which does a row-wise concatenation of the results of a function mapped to a list.

If you've found a solution to your problem that works for you (or almost works), feel free to add your own answer to your question. It's easier to work on and improve answers when it's done that way.

Thanks for all the responses, however, I got a relatable outcome from (https://unix.stackexchange.com/questions/113898/how-to-merge-two-files-based-on-the-matching-of-two-columns) by using
$ awk 'NR==FNR {h[$2] = $6; next} {print $1,$2,$3,h[$2]}' file_1 ref_f > out_1
For file_1 (NR==FNR is only true for the first file argument). Save column 6 in hash-array using column 2 as key: h[$2] = $6. Then for ref_f and output all three columns $1,$2,$3, appending the corresponding saved column from hash-array h[$2].
However, it does not print col4 and col5 from file_1 and must be repeated for each subsequent file. I hope that the answer could be improved further.


text:
What trajectory analysis method allows to set the form of the trajectory?
I have single cell data from two samples : normal and pathological and
I would like to track the progression of the cells (and find the genes that drive it) from normal to pathological.
Upon first glance, it seems most reasonable to define the trajectory in advance in such a way that it would have a starting point, an end point, and a straight route between these two.
Is it possible, and using what package if so?
P.S. In contrast to the question above - I started with the destiny Bioconductor package. However, it always yields three tips for the diffusion pseudotime, and I am not sure how to determine the starting and the end tip of the trajectory...

slingshot (Bioconductor) accepts a priori starts and ends. For differential expression along the trajectory you may check tradeSeq.