text:
What versions of pip and python are appropriate for nsdpy?
I am currently trying to download all the requirements to make a COInr database, following this workflow; https://mkcoinr.readthedocs.io/en/latest/content/overview.html.
It requires a couple different third party programs to be able to run, one of them being nsdpy, which it says all it requires is cloning the repo from Github and then you're supposed to simply be able to use 'pip' to install nsdpy.
However I continued to receive this error;
""DEPRECATION: Python 2.7 will reach the end of its life on January 1st, 2020. Please upgrade your Python as Python 2.7 won't be maintained after that date. A future version of pip will drop support for Python 2.7. More details about Python 2 support in pip, can be found at https://pip.pypa.io/en/latest/development/release-process/#python-2-support
ERROR: Could not find a version that satisfies the requirement nsdpy (from versions: none)
ERROR: No matching distribution found for nsdpy""
Which makes it seem that my issue is with the version of our Python (I am using a linux system on a Mac). The version of Python that is says we should be using is 2.7, the current version I have downloaded is 3.7.7. But then I was reading through the nsdpy docs (https://nsdpy.readthedocs.io/en/latest/main.html#overview) which says that it requires a Python version of 3.8+. So I am a little confused, which version of Python to I truly need to be able to install nsdpy?
Is there a more up to date version of nsdpy?
If anybody else has had this problem and/or has a solution, please let me know!

To add comment that was lost in migration: I strongly suspect that the version of pip is at fault, as it seems like the issue is that python2.7 is being run, not that it is required.


text:
Problem with Calculating Polygenic Risk Score (PRS) in Plink
I've been struggling with constructing a PRS score. In the pics you can see how my summary stats and my target data look like.
Target Data Summary Stats
I've tried the following code which should be pretty straightforward:
plink \
  --bfile target \
  --score summary_stats.txt 1 4 6 header \
  --out final

When I run it, I get this message:
10499994 variants and 1198 people pass filters and QC. Note: No phenotypes present.

The data is binary and the fam file really doesn't have phenotypes (it's -9). Is there any way to still run a PRS?
I'm unsure what I'm doing wrong and would really appreciate some insight! Also, the data has been harmonized already so it should technically work!

Please don't put images of text into posts. Copy-paste the actual text and add it as a code block. If formatting is difficult, someone else can edit the post to fix that up.

I'm a little confused about what you mean by there being no phenotype. To build the score you need a phenotype, although to apply it you might not. (Usually PRS is build with a base dataset and then applied using a target dataset). Here's a good tutorial about doing PRS in plink that might help you:
https://choishingwan.github.io/PRS-Tutorial/plink/
Also this paper is a really good paper that helped me better understand how building PRS works:
https://www.nature.com/articles/s41596-020-0353-1
Hope this helps!

source: Stack Exchange

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

sort: votes


text:
Best Library for 3D Protein Structure Visualization and Annotation
I have a background in Computer Science and have recently gotten very interested in the field of Protein Engineering. I am trying to build a web app which displays the 3D structure of a protein when I upload a PDB file. I also want to annotate subsequences with unique labels.
As I am very new to the field I am confused about which library to use. Can I use BioPython for everything? Or should I use any other library for 3D structure visualization like 3Dmol.js, NGL Viewer etc?
Any comments/suggestions/feedback would be really helpful guys as I am a beginner in the field of BioML.

BioPython can't do structure visualization. NGL and LiteMol are obsoleted by Mol* – from the same developers. There is venerable J(S)mol and its wiki lists a few other WebGL viewers. The newest viewer is Moorhen, but it's more for structural biologists.

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.

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


text:
How to annotate a bacterial genome automatically?
We recently got some shotgun sequencing results from a soil bacterium. And we have obtained some contigs which have some genes of interest.
Is there a way to automatically annotate the whole sequence (one single contig of around 25 kb size) and mark the genes on this stretch of DNA? I have never done this before so any help would be appreciated. The ultimate goal would be to know what genes are on this stretch of DNA and see if they work together in a pathway. But the first problem is:
How to annotate a stretch of DNA with many genes on it that would code for a protein/enzyme? so that I get a map of genes on this stretch of DNA.
Thanks.
Edit:
This is a soil isolate and we don't have a full genome sequence. Presumably, it is a Rhodococcus species, which might be similar. From the one gene we know, there are no exact homologs in the database so most probably it is a new strain. The strain works symbiotically with another strain which we are also working on and later I want to predict what proteins maybe involved in the catalysis of a specific compound. It is not a one step process so possibly multiple enzymes/proteins are involved.
I tried UGENE and UGENE can identify multiple ORFs in this particular stretch, but there are hundreds of ORFS (on both strands) and I find it tedious to blast on NCBI one by one.
I didn't want too many questions in one post, so at the moment just focus on the identification of all the putative genes on this one contig would be good enough to start.
Thanks again.

There are a number of annotation pipelines that you could use to get fairly standardized annotation - Prokka has been around for a while (paper published in 2014) but is still used frequently. There's also Bakta, which may be better at catching short ORFs and providing more detailed functional annotation, but requires a fairly large database.

After I know about Prokka, I found it as a service on Galaxy webservice here. https://usegalaxy.org/. Got the annotations done in no time.


text:
java.lang.ArrayIndexOutOfBoundsException: Index 86 out of bounds for length 86
I use the next versions:
gatk   4.4.0.0
minimap2   2.26-r1175
samtools   1.19

I have a fastq file. And try to implement gatk MergeBamAlignment to generate a new bam-file:
java -jar gatk.jar FastqToSam --FASTQ .fastq --OUTPUT u.bam
--SAMPLE_NAME sample1 -RG readgroup -PL ILLUMINA
minimap2 -a hg38.fa .fastq -o A.sam
java -jar gatk.jar SortSam -I A.sam -O A.sort.bam -SO queryname
java -jar gatk.jar SortSam -I u.bam -O u.sort.bam -SO queryname
java -jar gatk.jar MergeBamAlignment -O m.bam -R hg38.fa -UNMAPPED_BAM u.bam -ALIGNED_BAM A.sort.bam

But I receive error on the output of MergeBamAlignment:
INFO    2023-12-21 19:20:41     SamAlignmentMerger      Processing SAM file(s): [./myprogs/data/NR.Aligned.sort.bam]

[Thu Dec 21 19:20:42 2023] picard.sam.MergeBamAlignment done.

Elapsed time: 0.02 minutes.

Runtime.totalMemory()=486539264

To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp

java.lang.ArrayIndexOutOfBoundsException: Index 86 out of bounds for length 86
at htsjdk.samtools.util.SequenceUtil.calculateMdAndNmTags(SequenceUtil.java:990)
at picard.sam.AbstractAlignmentMerger.fixNmMdAndUq(AbstractAlignmentMerger.java:611)
at picard.sam.AbstractAlignmentMerger.mergeAlignment(AbstractAlignmentMerger.java:586)
at picard.sam.SamAlignmentMerger.mergeAlignment(SamAlignmentMerger.java:181)
at picard.sam.MergeBamAlignment.doWork(MergeBamAlignment.java:368)
at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:289)
at org.broadinstitute.hellbender.cmdline.PicardCommandLineProgramExecutor.instanceMain(PicardCommandLineProgramExecutor.java:37)
at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:160)
at org.broadinstitute.hellbender.Main.mainEntry(Main.java:203)
at org.broadinstitute.hellbender.Main.main(Main.java:289)

I have diagnosed u.sort.bam and A.sam.bam with ValidateSamFile:
java -jar gatk.jar ValidateSamFile -I u.sort.bam -M SUMMARY -R hg38.fa
No errors found
[Thu Dec 21 00:32:34 2023] picard.sam.ValidateSamFile done. Elapsed time: 0.00 minutes.
Runtime.totalMemory()=234881024

Tool returned:
0

java -jar gatk.jar ValidateSamFile -I A.sort.bam -M SUMMARY -R hg38.fa
ERROR    2023-12-21 01:17:28    ValidateSamFile    Requesting earlier reference sequence: 23 < 49
[Thu Dec 21 01:17:28 2023] picard.sam.ValidateSamFile done. Elapsed time: 0.18 minutes.
Runtime.totalMemory()=2113929216
To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp

Tool returned:
-1

I have done the same actions only with --MD in minimap2. But the same errors occur. I would like to find out how to fix this problem. I use the latest versions of the programs. I have downloaded hg38.fa each one time, but it didn't change the situation. I acts in according with gatk's article 'Data pre-processing for variant discovery' : https://gatk.broadinstitute.org/hc/en-us/articles/360035535912-Data-pre-processing-for-variant-discovery

I don't think that you need to resort as minimap should output the reads in the same order that they were ingested

However, the fact that the output of gatk SortSam isn't valid is disconcerting. Can you see if it's a small part of your fasta file that is the problem? perhaps chop it in half and see if you can reproduce it on a smaller fastq?


text:
Can I interpret 15x coverage area as a part of one strand is missing?
This is a BAM file from a 30x coverage Nebula Genomics test opened with IGV.
I found an area (1:30169255-30169600) where the coverage is less than 15x. Can I interpret it as a part of one strand of Chromosome 1 is missing?

UPDATE:
Seeing a different question where similar vocabulary is used: it would not be one strand that would be lost. single stranded DNA is not very happy in the nucleus, and would either lead to something seriously bad or be gap-filled to make it dsDNA again.
It would instead be that one of two chromosomes has a deletion at that position, which would be a classic CNV.
ORIGINAL ANSWER:
No.
There are a variety of copy number variant (CNV) calling tools that put a lot of effort into statistically robust ways of asking this question, all of which are a lot more complex than just looking at raw coverage. I'd recommend looking more into them.
In terms of explanations for lower coverage in some regions that aren't CNVs: sequenceability of DNA is not even, replication of DNA is not even, mappability of reads is not even, etc. You need to reject those hypotheses before you go to variants.


text:
Classical alignment with HISAT2
This question was also asked on Biostars
I have this script
#!/bin/bash

# Defining the number of threads to be used
threads=20

# Alignment of sequences to the reference genome using 'HISAT2'
echo -e "\nNow Running Classical Alignment-based with HISAT2...\n"

mkdir -p hisat2

# Download the human reference genome
wget -c ftp://ftp.ensembl.org/pub/release-100/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz -O hisat2/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
gunzip hisat2/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz

# Download the human annotation file
wget -c ftp://ftp.ensembl.org/pub/release-100/gtf/homo_sapiens/Homo_sapiens.GRCh38.100.gtf.gz -O hisat2/Homo_sapiens.GRCh38.100.gtf.gz
gunzip hisat2/Homo_sapiens.GRCh38.100.gtf.gz

# HISAT2 indexing to build the reference genome index
hisat2-build -p $threads hisat2/Homo_sapiens.GRCh38.dna.primary_assembly.fa hisat2/Homo_sapiens.GRCh38v3_hisat2.idx

# HISAT2 Alignment
for sample in $(cat sample_ids.txt)
do
    hisat2 -p $threads -x hisat2/Homo_sapiens.GRCh38v3_hisat2.idx \
        -1 trimmed_data/${sample}_1_val_1.fq -2 trimmed_data/${sample}_2_val_2.fq \
        -S hisat2/${sample}_hisat2.sam

    # Compress the SAM files to binary format, sort them, and index them.
    samtools view -Sb hisat2/${sample}_hisat2.sam | samtools sort > hisat2/${sample}_hisat2_sorted.bam
    samtools index hisat2/${sample}_hisat2_sorted.bam
rm hisat2/${sample}_hisat2.sam  # Remove the .sam files for storage purposes
done

# Quantifying Aligned reads using the Subread package's script 'featureCounts'
for sample in $(cat sample_ids.txt)
do
featureCounts -T $threads -a hisat2/Homo_sapiens.GRCh38.100.gtf \
        -o hisat2/${sample}_hisat2_counts.txt \
    hisat2/${sample}_hisat2_sorted.bam
done

# Generating analysis report with MultiQC
echo "Running MultiQC on HISAT2 alignment results ........."
multiqc -o multiQCreports hisat2
echo "MultiQC on HISAT2 alignment results completed successfully."

# Display message to show that the script is complete
echo "Script completed!"

However, I have receieved an error below
Load annotation file Homo_sapiens.GRCh38.100.gtf ...                       ||
||    Features : 1378275                                                      ||
||    Meta-features : 60683                                                   ||
||    Chromosomes/contigs : 47                                                ||
||                                                                            ||
|| Process BAM file SRR19170514_hisat2_sorted.bam...                          ||
ERROR: Paired-end reads were detected in single-end read library : hisat2/SRR19170514_hisat2_sorted.bam
ERROR: Paired-end reads were detected in single-end read library : hisat2/SRR19170514_hisat2_sorted.bam

Question: What can I do to solve this?

I believe the error is from featureCounts and the solution would be to add the -p flag to the featureCounts call, stating that you are providing paired-end data.
EDIT: You might want to check this r/bioinformatics (reddit) thread out.

Well, thanks. It worked when I added the -p flag to the featureCounts call

@EmmanuelAroma Accept haci's answer to mark the post as answered.

thank you so much comrades, I've got a wonderful answer


text:
Single Cell RNASeq Knee Plot - Collapsed vs Non-collapsed Reads
I am trying to adapt combinatorial index reads with some modifications for long read sequencing; it's a custom approach, not a commercial kit.
I have a question about using Knee Plot for identifying Barcode white list for non-10X data. While generating the Knee Plot barcode which barcode frequencies should be used: without collapsing PCR replicates to after collapsing PCR replicates?
I have used UMI-Tools (https://umi-tools.readthedocs.io/en/latest/QUICK_START.html) and Pacbio groupdedup based approach. For the knee plot generation I used both collapsed and non-collapsed reads. Since its long reads and coverage is not great; bumps are not the perfect one but both inputs generate it. I am not sure which input I should consider to remove the backgrounds.
Thanks
EDIT: Can the down-voters please leave a note why they think this not a good query to ask?

You're providing information as a trickle instead of providing all relevant information including what you've tried right in the post. One needs to read all comments to understand your position. Hence the downvote.

If you have UMIs, then cell/UMI/mapped transcript triplets should be collapsed into unique triplets.

I was looking at UMI Tools white list function; somehow it looks to me it consider all the reads for knee plot.


text:
Gene list from eclip-seq .bed files
I have some eclipseq data, three IP and three input datasets, and have the alignments for all three. I used Bedtools "Intersect Interval" with the genome gtf file (hg19) I used to align them with, and that creates a bed file with gene names for each peak. I then wrote an R script to sum the clipper score for each gene within a sample, then do a simple student's t test across groups. I now have a gene list and p values for differential expression, but not sure if this is correct.
I've tried looking at a lot of sources but cant really find any canonical way to generate a gene list from bed files.
I've also tried IDR and it doesn't work on the clipper output, unfortunately.
Thanks!
Edit:
eCLIP-sequencing, it's essentially a variant of iCLIP-sequencing and data analysis methods are transferrable, to my knowledge.
I have alignments (bam files) made with Star aligner against Hg19. After deduplication and sorting, I called peaks using Clipper (a peak-caller made by the creators of eCLIP). It outputs bed files with scores for each sample.
My question is how exactly do I turn this bed file with scores into a gene list? Is my method above correct?

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

What is "eclipseq data"? Can you please provide some example input files? Anything you can add to clarify the most important problem you're having will be useful for people who are answering this.

Added more information. Thanks!

As someone who works with eCLIP data from ENCODE, I think I can help out.
You need to think of this more similar to ATAC or RNA-Seq, whereby you actually get a counts matrix, normalize that, and then apply math.
Another good idea would be to look through this eCLIP analysis document provided by Gene Yeo of ENCODE for more specific ideas.
https://www.encodeproject.org/documents/1f171ac6-a36a-41ac-b632-741aeb47aad2/@@download/attachment/eCLIP_analysisSOP_v2.3.pdf

As others have said it's not very clear what you're asking here, but I'm assuming you want to know how to determine relevant genes that were enriched in your IP.
As per the group behind eCLIP, Van Nostrand et al. 2020
: "standard peak analysis used the set of peaks identified as irreproducible discovery rate (IDR) reproducible and meeting fold-enrichment (≥ 8-fold) and significance (p value ≤10−3) in immunoprecipitation versus paired size-matched input."
If those cutoffs are too stringent you could change them to fit the context of your project. And what do you mean by differential expression? It sounds like you're addressing binding.


text:
Seeking Advice on Acquiring Clean Experimental data comprises of proteins for FDR Estimation
I am a bioinformatics student doing my research on top-down mass spectrometry-based proteomics and seeking some guidance on acquiring very clean experimental data with known true positives and false positives. I want to assess my proposed FDR estimation method under controlled conditions and would appreciate your insights on obtaining high-quality experimental data.
My question is

Any tips on acquiring clean experimental data with known true positives?
Any suggestion on sample types or protocols you recommend for ensuring a reliable dataset for FDR estimation?
How can one introduce controlled errors or variations, like false positives, while maintaining the overall data quality? If anyone has any idea or experience or any paper suggestion in which they used a similar approach please share thanks.

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


text:
probability of finding a 5 amino acids in a row within a proteome
How to calculate the probability of finding two proteins that share a 5 amino acid long motif from a proteome of around 1067 proteins that have an average length of 65 residues.
The probability of a specific sequence that is completely the same is 1067/20^65

Its just significant at 2%

0.055 . 65 . 1067 = 0.021 (its actually 0.02 - see note)

Where critical probability is 0.05, i.e. 5%
There are 20 aa so 1/20 probability of random occurrence. The oligopeptide is 5 amino acids long, hence to the power 5. Obviously there are 1067 proteins and the probability is per protein, thus if these were independent its a straight multiplication because that is how independent events are measured.
Note The 65 is really 60 because it's a sliding window moving 1 amino acid at a time and the proteins are only 65 in length. The actual value is 60 rather than 65 because at the 3' terminal of the sequence the sliding window can't capture the remaining 4 amino acids, because once the oligopeptide < 5 the calculation is no longer valid. I used 65 so you could see where the value is originating.
Obviously the proteins in question are not independent so you need to state the assumption, i.e. assuming independence.

There are 22 amino acids, not 20. Almost all animals (the only exceptions I know are a few insect species) have 21, because of selenocysteine, the 21st amino acid. I know less about pyrrolysine, the 22nd, but I believe it is found in some bacteria and archaea only. Given how few occurreneces of this amino acid most species have, I doubt it will change the analysis much, but there are 21 amino acids in most non-plant eukaryotes.

Thanks @terdon, now I know.  Selenocysteine is interesting because of the importance of cysteine disulphide bonds. I assume selenocysteine doesn't form them. Anyway, the probability of random occurrence is 1/21 or 1/22 depending.

I... don't remember. Wow. I did my PhD on these things but it's been a while and I always focused on detecting them rather than on their biochemistry, but still! Anyway, Sec has selenium where Cys has sulfur, so presumably it cannot make disulphide bonds by definition, but I do seem to recall that the selenium allows for very similar bonding. We tended to treat Sec as a more reactive version of Cys, but my ignorance of biochemistry is deep so I don't really know.


text:
Process gets killed before completing
I have a large tsv file (~15 GB) that is a kmer count matrix and I am trying to calculate jaccard distance of the counts but every time I try to run the following code the process ends before it's done with a killed error message in the STDOUT:
jaccard = pd.DataFrame(squareform(pdist(kmer_matrix, 'jaccard')))
jaccard.to_csv(matrix_path+kmer+"_"+method+"_jacDistance_noFilter.tsv", sep='\t')

This process works really well with smaller tsv files so I know the code works but I wanted to know if anyone knew what would help run the process quicker/completely. I know the problem must be with my RAM/CPU but is there a way to maybe break up the matrix to calculate everything in batches to make it easier to work with?

You've overloaded the memory. One route is to deploy "chunking". pandas will do this. This is complicated to implement code wise however.
My suggested solution (because it simpler) would be to do this via pandas2 which uses Apache Arrow via PyArrow. You need to specify which method is being used at the onset of declaring the dataframe. Apache Arrow is for lifting large amounts of data. If this works there's no complicated coding you just import pandas2 as pd and in the dataframe construction describe the precise heavy lift method (I think there's an additional approach besides PyArrow).

Thank you! I was able to use a chunking method to calculate all the distances! Now the issue I'm having is I can't open the file to see the distances. Normally, I'd convert the .tsv file to a .meg file by going into the file and adding the proper headers to visualize the tree with MEGAX. But since I can't open the file I can't normally convert it. Any tips?

Hi @rimo, great! Can you put this as a separate question please? It should be easily solvable. If you do I'll answer it on Saturday, because I've work commitments.

Thanks @M__!! I appreciate your help! I just posted another question. Take a look when you can!


text:
Fastq file with very high per base sequence quality and no box-whisker plot
I have my data from NCBI. I am having difficulty making sense of and correlating my fastq file read 1 and read 2 and the fastQC report. As shown below, I have two paired-end read results each displaying 127 bp in length; however, my results on the fastQC report of per base sequence quality show 91 and it does not show the interquartile range (box-whisker plot). What could be the problem?
Below is an excerpt display of my fastq file and then the fastQC report:
@SRR20708369.1 NS500749:581:HGV7VBGXB:1:11101:17225:1048 length=127
AGCTATCACAGCANGGTACTTCCCAGCGTTACGGGANTGGCCGTACGCCTAACCGCTAACATTACTGCAGGCAACCTACTCATGCACCTAATTGGAAGCGCCACCCTAGCAATATCCACCATTAACC
+
A6AAAE6EA/AAA#EAEAEEEAEEEE/EAEEEEA<E#AAAAEEEE6E<E//EE6EAEEAEAEE/AAA6E//E/EE/A<AEEEE/AEEE//<EEA/EE/E/6<<EE<EEE/<EEEEE//<AE/AAEA/

@SRR20708369.2 NS500749:581:HGV7VBGXB:1:11101:22488:1048 length=127
AGCTATCACGCCANAAGGCTTCCGAGTTAGTATAGGNGCAAAGTGAATACATCTCTGTAACCAGACTCAGATCAAGAAATGATACATTCTGGCCAGGTGTGGTGGCTCATGCCTGTAATCTCAGCAC
+
AAAAAE/EAAAAA#EEEEEEEEEEEEEEEEEEEEE6#AA/AEA/EEEEEEEEAEEEEAEA/EAE/AEAEEA6</6E///EE//6A//AEAEE<E/EE<E<EEEEEEEEEAEAEAAE///EEE/EE6/

:

Cross-posted and answered on biostars: https://www.biostars.org/p/9583849/#9583850

https://trace.ncbi.nlm.nih.gov/Traces/?view=run_browser&acc=SRR20708369&display=metadata is a 10x single-cell sample.
So FastQC plot you are showing us appears to be the R2 read from this sample which is 91 bp (split files, where as the example read you show above seems to be Illumina Index + R1 + R2). This is normal. Illumina index (8 bp), Read 1 is 28 bp (Cellbarcodes + UMI) and R2 is RNA (91 bp) = 127 bp.


text:
LD local clump-GRCh38 bfile
I want to do LD local clumping.
The tutorial (https://mrcieu.github.io/ieugwasr/articles/local_ld.html) teaches how to, but the downloaded bfile is GRch37.
I followed this question to download GRch38.
LD clump GRCh38 GWAS results
I downloaded and it is "all_hg38.pgen". How can I transform it into bfile? Thank you.

I have the same problem when facing the analysis of eQTL data for clumping and I do not find any answer to solve it. So I try to find a method to transform 1000 Genomes phase 3 ailgned by hg38 data which is pfile format (pgen, pvar, psam). I found a method at: https://cran.r-project.org/web/packages/snpsettest/vignettes/reference_1000Genomes.html. This is a method about how to extract the EUR population and transfrom to bfile format using PLINK 2.0. After this transfrom we can use PLINK 1.9 to clump the GWAS data. I think this is the method which we want to use and I will try it. Hope this can help you.


text:
why would a pdb ID generate 2 templates that are identical
I used the sequence of a protein kinase to generate templates for another and build a homology model using these templates using swiss-model online tool. there are 2 templates that were generated using the same pdb ID and I cannot find a reason as to why this happened, they are identical in all aspects and the homology models generated were also identical with no variance

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.

please write down in your question all the passages you did starting from here: https://swissmodel.expasy.org/interactive

have you read the documentation ? https://swissmodel.expasy.org/docs/help apparently templated are not generated but found given your input , what does 'The Template Results page' serves to you ??

What is the sequence that you used? What templates were produced?


text:
problem with feature counts of RNA bulk seq paired data- in Rsubread
I did manage to bild index file (from NCBI transcriptome) and perform alingment in Rsubread in R with my fq files. I did get BAM files as a result and no error. However i am having trouble with next step as featurecount function does not produce any result( 0%). I have tried builtin annotation (hg38)
i am wondering if BAM files are ok..
this is the code i have used
library(Rsubread)
index_path <- "D:/FASTQ_Raw/hum_transcr"
all_files <- list.files(path = "D:/FASTQ_Raw", pattern = ".fq$", full.names = TRUE)
common_prefixes <- unique(gsub("_[12]\\.fq$", "", basename(all_files)))

output_directory <- "aligned_reads/"

dir.create(output_directory, showWarnings = FALSE)

for (prefix in common_prefixes) {

  read1_file <- paste0("D:/", prefix, "_1.fq")
  read2_file <- paste0("D:/", prefix, "_2.fq")

  output_file <- paste0(output_directory, prefix, "_aligned.bam")

  align(index = index_path,
        readfile1 = read1_file,
        readfile2 = read2_file,
        output_file = output_file,
        nthreads = 4)
}

so the results/stats of individual BAMs appear ok still when i try featurecount function on BAM files i get zero....i was also using built in annotation index hg38 for attempt of counting but it does not work. I had issues with getting  gtf file from NCBI. What do you recommend ?

How did you run featureCounts() and what was the output? Are you sure your reference and gtf file are compatible?

Well i have created Index file in first step from .fa transcriptome file from NCBI ( for gh38.14). Howeverr there i had problem of opening gtf file. I did download one from ensembl but i am not sure if they are 100 % compatible ( they should be i guess). Because of this i have tried Builtinannotation function from Rsubread but this yielded 0% alignment. I was hoping to do everything in R because i am beginer with seq analysis and with R I have at least some experience with different packages.

NCBI and Ensembl use a different format for chromosome names. That's probably why you didn't get any results. The fasta and gtf must be from the same dataset.

Many Thanks!! I will retry with NCBI gtf file. Alignment was quite reasonable (over 90%) so it should work...

I don't think FeatureCounts is meant to work on an alignment to the transcriptome.  It's not at all clear to me that the subread aligner, or any aligner, is meant to be used on a transcriptome.
If you have the transcriptome, I'd just use pseudoaligners like Kallisto or Salmon.

well i think it does, at least according to the user manual. And feature counting is mandatory step after alignment to the index. Not that i am an expert .IS Kallisto Windows compatible? I am not familiar with Linux...

Aligners can definitely be used on a transcriptome. That's done all the time on non-model organisms for which full genome is not available. I don't think the Rsubread aligner is different in that sense.

But they aren't going to handle ambiguous alignments properly, like a psudoaligner would.


text:
conducting a ttest on a board of 1000 genes in 50 peope divided into to groups(A and B). pvalues are different each time
I am struggling in conducting at test on a board of 1000 genes in 50 individuals divided into to groups (A and B). I have a dataset of 1000 genes of two groups of people (total individuals=50)
group A = 30 (healthy) individuals & group B = 20 sick) individuals
I created a matrix where rows are genes and ind are columns.
now I need to perform a t-test for all the genes and find the p values of my alternative hypothesis (not the null one): i)A<B, A>B and A=!B(not equal). threshold is 0.05 after FDR (which I must include)
I am struggling with the code since t-test differ on sensitivity and each of my approach returns a different result. The dummy data is as follows:
import pandas as pd
import numpy as np

# Generate data for group A (healthy individuals) & group B (sick individuals)
A = np.random.rand(30, 1000)
B = np.random.rand(20, 1000)
# Create dataframe for group A & B with genes as columns
Adf = pd.DataFrame(A, columns=[f'Gene_{i}' for i in range(1, num_genes + 1)])
Bdf = pd.DataFrame(B, columns=[f'Gene_{i}' for i in range(1, num_genes + 1)])

# Add a column indicating the group
Adf['Group'] = 'A'
Bdf['Group'] = 'B'

df = pd.concat([Adf, Bdf], ignore_index=True)
df.head()

Any ideas?

can you please post the code that you've tried so far, along with a minimal example of your data? what are the approaches you've tried and why is it unexpected that different approaches yield different answers?

sure. thanks for the response. so basically i have a dataset that i turned into a python matrix (with pandas DataFrame) consisting of 50 individuals, 30 individuals in group A (the first 30 columns) and 20 individuals in group B. the individuals in group A are healthy and the individuals in group B are sick. i applied the t-test for all genes but i have to find the pvalues for Hypothesis 1 (alternative) if i) A < B, ii) A > B, iii) A != B (not equal); threshold is 0.05 after FDR
(code at my next comment)

group_A = dataset.iloc[:, :30]  # first 30 columns for gr A (up to 30, including 30)
group_B = dataset.iloc[:, 30:]  # next 20 columns for gr B (next 20, from 31 to 50)
p_values_less = []
p_values_greater = []
p_values_not_equal = []
t-tests for each gene
for gene in range(dataset.shape[0]):  #access to all genes
t_stat, p_val_not_equal = stats.ttest_ind(group_A.iloc[gene], group_B.iloc[gene])
p_values_not_equal.append(p_val_not_equal)
p_values_less.append(p_val_not_equal / 2 if t_stat < 0 else 1)
p_values_greater.append(p_val_not_equal / 2 if t_stat > 0 else 1)

p_values_corrected_less = smm.multipletests(p_values_less, alpha=0.05, method='fdr_bh')[1]
p_values_corrected_greater = smm.multipletests(p_values_greater, alpha=0.05, method='fdr_bh')[1]
p_values_corrected_not_equal = smm.multipletests(p_values_not_equal, alpha=0.05, method='fdr_bh')[1]
results = pd.DataFrame({
'Gene': dataset.index,
'p_value_less': p_values_corrected_less,
'p_value_greater': p_values_corrected_greater,
'p_value_not_equal': p_values_corrected_not_equal
})

significant_statistic_results = results[(results['p_value_less'] < 0.05) |
(results['p_value_greater'] < 0.05) |
(results['p_value_not_equal'] < 0.05)]
print(significant_statistic_results)

it messed my lines and Tabs but thats my code

Please [edit] your question and add the code there. You can use the {} button to format it as code.

group_A = dataset.iloc[:, :30]  # first 30 columns for group A (up to 30, including 30) group_B = dataset.iloc[:, 30:]  # next 20 columns for group B (next 20, from 31 to 50) #hypothesis:i)Α < Β, ii) Α>Β, iii) A != B p_values_less = [] p_values_greater = [] p_values_not_equal = [] for gene in range(dataset.shape[0]):     t_stat, p_val_not_equal = stats.ttest_ind(group_A.iloc[gene], group_B.iloc[gene])     p_values_not_equal.append(p_val_not_equal)     p_values_less.append(p_val_not_equal / 2 if t_stat < 0 else 1)     p_values_greater.append(p_val_not_equal / 2 if t_stat > 0 else 1)

`# Apply FDR correction
the p values for each hypothesis (H1) are: p_values_less, p_values_greater, p_values_not_equal
#a = 0,05
#FDR=b/h
p_values_corrected_less = smm.multipletests(p_values_less, alpha=0.05, method='fdr_bh')[1]
p_values_corrected_greater = smm.multipletests(p_values_greater, alpha=0.05, method='fdr_bh')[1]
p_values_corrected_not_equal = smm.multipletests(p_values_not_equal, alpha=0.05, method='fdr_bh')[1]`

`results = pd.DataFrame({
'Gene': dataset.index,
'p_value_less': p_values_corrected_less,
'p_value_greater': p_values_corrected_greater,
'p_value_not_equal': p_values_corrected_not_equal
})
significant_statistic_results = results[(results['p_value_less'] < 0.05) |
(results['p_value_greater'] < 0.05) |
(results['p_value_not_equal'] < 0.05)]
print(significant_statistic_results)`

Could it be that you are getting different results because you are not setting a seed on your random number generator?
The recommended way to make sure you get the same results when using numpy to generate dummy data is:
import numpy as np
# make random number generator with seed
rng = np.random.default_rng(seed=42)
rng.random([4,4])
#Out:
#array([[0.77395605, 0.43887844, 0.85859792, 0.69736803],
#       [0.09417735, 0.97562235, 0.7611397 , 0.78606431],
#       [0.12811363, 0.45038594, 0.37079802, 0.92676499],
#       [0.64386512, 0.82276161, 0.4434142 , 0.22723872]])

```


text:
Fetching Protein Sequences through UniPROT gives ConnectionError
I have 193000 protein interactions in CSV named all_proteininteractions.csv which have query protein name in the 2nd column and partner protein in the 3rd like this:
Query_ENSP,Query_Name,Partner_Name,Combined_Score
9606.ENSP00000000233,ARF5,RAC2,0.4
9606.ENSP00000000233,ARF5,RAC2,0.4
9606.ENSP00000000233,ARF5,PPP2CB,0.408
9606.ENSP00000000233,ARF5,PPP2CB,0.408
9606.ENSP00000000233,ARF5,R3HDML,0.413
9606.ENSP00000000233,ARF5,R3HDML,0.413
9606.ENSP00000000233,ARF5,PLEK,0.422
9606.ENSP00000000233,ARF5,PLEK,0.422
9606.ENSP00000000233,ARF5,RAB11FIP5,0.501
9606.ENSP00000000233,ARF5,RAB11FIP5,0.501
9606.ENSP00000000233,ARF5,AP3B1,0.553
9606.ENSP00000000233,ARF5,AP3B1,0.553
9606.ENSP00000000233,ARF5,CYTH4,0.579
9606.ENSP00000000233,ARF5,CYTH4,0.579
9606.ENSP00000000233,ARF5,RABAC1,0.658
9606.ENSP00000000233,ARF5,RABAC1,0.658
9606.ENSP00000000233,ARF5,YKT6,0.672
9606.ENSP00000000233,ARF5,YKT6,0.672
9606.ENSP00000000233,ARF5,BET1,0.679

I am using this Python code to fetch sequences of each protein and store them in a CSV file:
import pandas as pd
import requests
import math

requests.packages.urllib3.disable_warnings()

# Function to fetch protein sequence from UniProt
def fetch_protein_sequence_from_uniprot(protein_name):
    #Declaring UniProt API
    uniport_api_url = f"https://rest.uniprot.org/uniprotkb/search?query=gene_exact:{protein_name}+AND+organism_id:9606&format=fasta&size=1"
    response = requests.get(uniport_api_url,verify=False)

    # Parse the response to extract sequence
    sequence = ""
    if response.ok:
        lines = response.text.split("\n")
        for line in lines:
            if not line.startswith(">"):  # Exclude header lines
                sequence += line
    return sequence

# Reading the input CSV file
input_csv = "all_protein_interactions.csv"
chunksize = 1000  # Batch size

# Create an iterator to read the CSV in chunks
csv_reader = pd.read_csv(input_csv, chunksize=chunksize)

# Initialize an empty list to store modified chunks
modified_chunks = []

# Process each chunk
for chunk in csv_reader:
    # Fetch protein sequences for query and partner names
    query_sequences = [fetch_protein_sequence_from_uniprot(name) for name in chunk.iloc[:, 1]]
    partner_sequences = [fetch_protein_sequence_from_uniprot(name) for name in chunk.iloc[:, 2]]

    # Add fetched sequences as new columns in the chunk
    chunk['Query_Sequence'] = query_sequences
    chunk['Partner_Sequence'] = partner_sequences

    # Append the modified chunk to the list
    modified_chunks.append(chunk)

# Concatenate all modified chunks into a single DataFrame
result_df = pd.concat(modified_chunks, ignore_index=True)

# Save the updated DataFrame back to CSV
result_df.to_csv("protein_interactionswithseqs.csv", index=False)

First, it was throwing SSLError, so I made verify=False. But now, after running the script for 3 hours it's showing following error:
ConnectionError: ('Connection aborted.', RemoteDisconnected('Remote end closed connection without response'))

How can I remove this connection aborted error to fetch 386000 sequences?

You also should add print statements in various places, to see where the problem is. Has the script managed to finish the for chunk loop? Where exactly is it failing? And isn't your loop trying to get the same protein multiple times? You seem to have the same name in the 2nd column for several lines, so why get that same protein multiple times?

Thanks. I am using  Windows 10 with Intel(R) Core(TM) i7-3770 CPU and 12 GB ram.

I modified the text as well. I have similar proteins but I want to fetch them all and store in 5th and 6th column

Your script works as expected on the example you gave. Running your script results in a protein_interactionswithseqs.csv with lines like:
$ head -n2 protein_interactionswithseqs.csv
Query_ENSP,Query_Name,Partner_Name,Combined_Score,Query_Sequence,Partner_Sequence
9606.ENSP00000000233,ARF5,RAC2,0.4,MGLTVSALFSRIFGKKQMRILMVGLDAAGKTTILYKLKLGEIVTTIPTIGFNVETVEYKNICFTVWDVGGQDKIRPLWRHYFQNTQGLIFVVDSNDRERVQESADELQKMLQEDELRDAVLLVFANKQDMPNAMPVSELTDKLGLQHLRSRTWYVQATCATQGTGLYDGLDWLSHELSKR,MQAIKCVVVGDGAVGKTCLLISYTTNAFPGEYIPTVPGRPCALLQDAEAGLSSDVEGACWCQTQEGRMTGEAPAPPWPALRGPHHS

However, this is horribly, horribly inefficient! You are fetching the same sequence multiple times. With your example data, the protein ARF5 is actually downloaded 19 times, meaning this will take 18 times more time than needed. This is very likely the cause for the timeout error you describe. Now, it might be that your sequences are too many to be stored in RAM (but I would still try it), so one way around this is to download the sequences first, and then write the mappings. Something like this:
import pandas as pd
import requests
import sys
requests.packages.urllib3.disable_warnings()

def write_protein_seq_file(protein_name, protein_seq):
    file_name = "seqs/" + protein_name
    seq_file = open(file_name, "w")
    data = f"{protein_seq}\n"
    seq_file.write(data)
    seq_file.close

def read_protein_seq_file(protein_name):
    seq_file = open("seqs/" + protein_name, "r")
    data = seq_file.read()
    seq_file.close
    return data.rstrip()

# Function to fetch protein sequence from UniProt
def fetch_protein_sequence_from_uniprot(protein_name):
    #Declaring UniProt API
    uniport_api_url = f"https://rest.uniprot.org/uniprotkb/search?query=gene_exact:{protein_name}+AND+organism_id:9606&format=fasta&size=1"
    response = requests.get(uniport_api_url,verify=False)

    # Parse the response to extract sequence
    sequence = ""
    if response.ok:
        lines = response.text.split("\n")
        for line in lines:
            if not line.startswith(">"):  # Exclude header lines
                sequence += line
    return sequence

# Reading the input CSV file
input_csv = "all_protein_interactions.csv"
chunksize = 1000  # Batch size

# Create an iterator to read the CSV in chunks
csv_reader = pd.read_csv(input_csv, chunksize=chunksize)

# Initialize an empty list to store modified chunks
modified_chunks = []

# Get all sequences
for chunk in csv_reader:
    for protein in set(list(chunk.iloc[:, 1]) + list(chunk.iloc[:, 2])):
        protein_seq = fetch_protein_sequence_from_uniprot(protein)
        write_protein_seq_file(protein,protein_seq)

print("Finished fetching sequences", file=sys.stderr)
## Now, process the file again and write the new file
csv_reader = pd.read_csv(input_csv, chunksize=chunksize)

for chunk in csv_reader:

    # Fetch protein sequences for query and partner names
    query_sequences = [read_protein_seq_file(name) for name in chunk.iloc[:, 1]]
    partner_sequences = [read_protein_seq_file(name) for name in chunk.iloc[:, 2]]

    # Add fetched sequences as new columns in the chunk
    chunk['Query_Sequence'] = query_sequences
    chunk['Partner_Sequence'] = partner_sequences

    # Append the modified chunk to the list
    modified_chunks.append(chunk)

# Concatenate all modified chunks into a single DataFrame
result_df = pd.concat(modified_chunks, ignore_index=True)

# Save the updated DataFrame back to CSV
result_df.to_csv("protein_interactionswithseqs.csv", index=False)

First, create a directory named seqs, then run the script above. It will process the csv file twice. Once to get all sequences mentioned in the file and save each one in a file called seqs/$sequence_name (where $sequence_name is replaced by the actual name of each protein. Then, it goes through the file a second time, but instead of the very time consuming download, it now fetches the sequence from the file stored on your computer. Once you are satisfied this works, you can delete the seqs directory.
Note that I don't have access to a Windows machine, so this was tested on Linux. Also, I am not a python coder, so I am sure this ins't very good python and there will be improvements to be made, but even on the toy example you gave, this speeds things up significantly:
$ time your_script.py

real    0m10.186s
user    0m0.998s
sys     0m0.076s

$ time my_script.py
Finished fetching sequences

real    0m3.180s
user    0m0.456s
sys     0m0.066s

As you can see above, even on the tiny example with only 19 data lines, your script took 10 seconds while mine only took 3. So try this approach (or, better, take my example and make it more pythonic, add some error checking, some try/except blocks, make sure your data are correct) and you might find the timeout error is no longer an issue. You should certainly find that the whole thing takes far less than three hours(!) to run.
First of all, that is a really, really inefficient way of doing it. Using your example data, you are getting the same protein, ARF5, 19 times. This makes your program 18 times slower than it needs to be. You should modify the script to download each protein once and then write them all out.

Your example script is not taking any file and where to download the all sequences

@AnasJamshed what do you mean? It is taking the exact same file as your version and in the exact same way. The sequences will be downloaded into the seqs directory which, as I explain in my answer, you need to create.


text:
Suggestion for way foward with upregulated genes
I am working with single cell RNA data- control and treated. I analyzed and got a list of upregulated genes. Upon doing pathway analysis all the genes were of apoptotic pathway. But these genes can signify different meaning to my interpretation, because, these genes were upregulated even in diseased condition. let me take an example, NEAT1 is one of the gene that is seen upregulated in the cancer condition and is associated with cancer progression. the treated sample which is a therapy also has neat1 upregulated compared to control which was untreated. Just because pathway analysis showed apoptosis, I cannot conclude that neat1 gene is aiding in killing the cancer cell due to the therapy right/ It can also be like due to continued upregulation its delaying cancer cell death due to therpay.
How to confirm computationally the above two statements and find the upregulated genes actual role in the therapy?
to be noted: there is no prior data or pattern available for the therapy based gene expression yet.


text:
Where can you download functional domain labels of the human genome?
I'm visualizing where mutations occur in a gene, and would like to add functional domains to the plots.
I have some trouble locating sources to download data on functional domains.
The dataset should contain:

Functional domain name
Start position in gene
End position in gene
Gene symbol name

Any tips kindly appreciated!

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

I would suggest looking into GFF3 annotations of the human genome.
If what you want is protein domain annotations, you may be more interested in the UniProt database which does seem to include genome annotation tracks for different protein functional regions.
I am unclear if uniprot directly has the gene names but this would be trivial to generate with a bedtools intersect operation.


text:
What is a common number of cells per cell type in single cell RNA Seq?
Good afternoon,
I just started working with Single cell RNA Seq and I am trying to understand what is an average / possible number of cells per sample.
I am looking at a dataset of 7 patients with kidney tumour samples (please find attached). Would these number of cells be possible? I basically want to double check whether I didnt do any mistakes in my pre-processing.
Thank you

It depends on many factors, both biological and technical, and it's impossible to answer in general.
The starting question is which technology was used, and how many cells were used in input. For example, in a SMART-Seq experiment with one cell per well (in a multi-well plate), you should know how many cells you have (one per well). For a 10x experiment, the person running the experiment should have measured the cell concentration in the suspension before loading the chip, and targeted a particular number of cells (typically 10,000). In a combinatorial barcoding approach (e.g. Parse), you'd also have a measurement beforehand so as to load a particular number of cells per well. So, you should already have an idea of the total number of cells.
Then, the first step of the analysis is typically some QC. In barcode-based approaches (including 10x), you'd make a knee-plot to estimate the actual number of sequenced cells in your sample.
In the clustering and annotation steps, you'd start distinguishing and annotating groups of cells, which hopefully do correspond to cell types. The algorithms are not perfect, so the clusters don't necessarily match cell types perfectly. For example, a rare cell type will easily be lumped together with another, similar, cell type. For a well-represented cell type, the algorithms have more chance to capture subtypes or cell states.
And of course there is the biology, and your knowledge of it. How many astrocytes are there in the kidney? I assume any number above 0 is worrisome. I'd expect the cell types in the data to roughly reflect the cell types in the tissue (but not exactly, there can be all kind of bias).
There is a question of depth: if you have few cells in a cluster, and few reads per cell, the gene expression data in that cluster is not very reliable. The actual number depends on the technology and read depth. I would say for 10x I would prefer at least 100 cells, and would not trust a cluster with <10 cells.
So for your question, you need to use your knowledge of the kidney and of the experimental design to decide. You can also compare to published scRNA-Seq experiments of the same tissue. You can also compare between patients to see if there is something troubling (for example, why does patient 1 have epithelial cells but no adipocytes, and opposite for patient 7? It might not be a problem, but only with context can you tell). And importantly, you do need to check that the annotation does match known gene expression patterns.

Dear Alexlok, Thank you so much for your detailed and very helpful answer, I will save it in my notes :). I have read a couple of times now the sentence "use your knowledge of biology" e.g. the found cell clusters. I definitely need to increase my cell biology knowlege here I think.... And yes also consider the experimental design more.


text:
What is the best way to acquire protein isoform sequence alignment?

[Update] Thanks @terdon. To clarify my question:
I have a bunch of protein isoforms sequences (produced by the transcripts in the figure) and I want to align protein products (e.g., proteins produced by the first 3 protein-coding transcripts against each other). See there are 3 blocks missing in 5' region in BRCA2-204 compared against BRCA2-206.
This case might be very easy with any MSA softwares such as clustal omega, but there are complicated cases which may trick software into aligning single amino acids instead of being splicing aware. Instead of manually adjusting the alignments, I want to know whether I can directly download, or "translate", this figure into alignment results (e.g., fasta format or clustal format)
This figure should be visible with https://www.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=ENSG00000139618;r=13:32315086-32400268
in the gene summary tab
[Original] Is there a way to directly download the sequence alignment result from Ensembl? If not, what will be the best way to align exons and introns without parameter-tuning

Please [edit] your question and clarify what you need. You mention protein isoforms in the title, but show transcript sequences in the body. Do you want the alignment of the proteins produced by each transcripts against each other? The alignment of each transcript against a reference genome? The alignment of the transcripts against each other? Ideally, give us an example we can look up, a link that we can click on to see what you are using (is that human BRCA2 or something else?) and then explain what you would like to get from that.

@terdon Thanks!

Interesting question (upvoted). This might be naive, but if these are transcripts surely they have already been spliced? That aside from eukaryotic genomics its easy, but must have been solved years ago, I suspect a de novo pipeline is simply "reinventing the wheel".

@M__♦ I am not sure I followed. You mean aligning these spliced exons is easy?

I've undeleted the question ... it's a perfectly good question - it has an answer by an expert which is upvoted. My justification is that I had trouble logging in over the last 24 hours - I had to escalate to SE login HQ to get in, hence I couldn't take part in the discussion before the question was deleted. @terdon is an expert on human genetics and I'll defer to their answer.

Please don't deface your question. Once a question has been posted and answered, you have given SE the right to use the content under CC-BY-SA. Even if your question is no longer useful to you, or if you have your answer, that doesn't mean that other users won't also find it useful so deleting useful content from the site isn't allowed.

Here is a manual (or semi-manual, at best) approach. If you need to do this for very many sequences it won't be practical, but if dealing with one or two genes, it might be enough. The idea is to:

Get the protein coding transcripts
Get the proteins they code for
Align

In the page you gave, there is a "Show transcript table" button near the top:

If you click on that, you get the list of transcripts and their corresponding proteins:

You can then click on the "Show/hide columns" and select "Translation ID" so you can get the Ensembl protein corresponding to each transcript:

If you now deselct everything and only leave "Translation ID" selected:

You are left with a nice list of Ensembl proteins. Copy/paste those into a text file:
ENSP00000369497.3
ENSP00000499438.2
ENSP00000505508.1
ENSP00000439902.1
ENSP00000514856.2
ENSP00000434898.2
ENSP00000506251.1
ENSP00000499570.2
ENSP00000433168.2
ENSP00000499256.2
ENSP00000514855.1

Now, go to BioMart and retrieve the sequences. Select "Ensembl Genes" as the database and then "Human Genes":

Then, click on "Filters" (1) and then "Gene" (2), select "Input external references ID list [Max 500 advised]" (3), set the ID to "Protein stable ID(s) with version" (4) and paste the protein list (5):

Next, go to "Attributes" (1) and select "Sequences" (2):

Now, finally, click on "Results" (1) and you will get the sequences in fasta format (2):

I'm afraid they are too long to fit in this answer, but I posted them here: https://pastebin.com/NZycAmgD.
You can now paste these sequences into your favorite aligner and get the alignment you were after. Here, I used the online clustalΩ (the image only shows the first few lines):

You can find the whole alignment in clustal format (text) here.

Thanks @terdon. This is exactly my question. In your clustal omega result screenshot (last figure of your answer) there are a bunch of misalignments for ENST00000700201 compared with others (starting from loc ~105 to the end). I assume the software is forcing different exons aligned together to reduce number of gaps. I mean those gaps should be inserted before loc ~105 on ENST00000700201, not after. But somehow the software prefer before (due to parameters, misalignments outperformed gaps). I can manually adjust it. Easy for one gene, but not for hundreds of genes. That's why I want to downloa

although ENST00000700201 is NMD, not a perfect example. But the question stays for other protein coding transcripts.

@JenMarylin I don't quite understand what you mean. What exons? We're aligning simple protein sequences here, the aligner has no knowledge of the exonic structure of the transcripts that were translated to make these proteins. And this is a multiple sequence alignment, it is positioning the gaps so as to maximize the score of the entire alignment. What is it you want to download? Are you looking for an alignment of the translated exons instead of the actual protein product?


text:
Approaches for filtering bacterial/fungal contaminant sequences from RADseq results
I'm working with RADseq data (288 compressed FASTQ files, *.fq.gz) from 3 plates of plant tissue. I've done demultiplexing and preliminary QC (using the Stacks process_radtags command), and before performing a reference alignment, I want make sure the sequences in each sample (roughly) correspond to the organism I expect them to be derived from, rather than potential contaminants (bacteria and/or fungi).
Does anyone have a straightforward, "canned" approach for doing this, either in R or a different environment? I believe the blastn algorithm can be used to query the sequences in these files against a database (like NCBI), and is usable on the command line. However, I've not used it before, and I'm not sure if it makes sense to query all the sequences in each file. Additionally, I don't know if the command can process this many files of this size (~600 MB per FASTQ file).
This post from ~4 years ago asks a similar question for a pair of human blood samples, but is concerned with bacterial/viral DNA. The answer linking kmerfinder might have promise, but I'd appreciate any input from folks familiar with this issue.
(Background: I've gotten my Master's degree in Biology and have worked with NGS data for ~4 years now. I'm performing analyses on a 32-core Linux server. I'm less skilled in Python, but am willing to use it if there are no other options.)

Do you have a high quality reference genome of the plant you are working with?

Yes, I do. I've used the same reference genome in other projects, with similarly related taxa.

Thanks for the question, I'm familiar with working on RAD format and particularly don't know whether you've generated mixed genome RAD data. If know this information and I've time I'll put an answer down using either Bash or Python. I don't use R unless I've no other option. Please post a note below the question if you do edit it to notify me.

I ultimately addressed this issue using the kraken2 software. I found this approach to rely on a good deal of my server's memory, but it was sufficient for doing a viral/bacterial classification of sequences.


text:
How to adjust cutoffs in Mutect2
I'm currently performing some initial analyses on targeted Illumina sequencing data generated with DNA from diseased tissue. My hypothesis is that one particular gene might contain somatic mutations in this disease. I'm currently running Mutect2 in tumor-only mode, using default settings.
As I have only found very few mutations (this result might well be plausible), I'd like to adjust the Mutect2 cutoffs to increase the sensitivity and see whether the result changes. For instance, I'd like to reduce VAF, depth of coverage and QUAL. Unfortunately, I can't work out which of the many arguments have to be modified and how: https://gatk.broadinstitute.org/hc/en-us/articles/360037593851-Mutect2
Can you help me out? Thanks!
gatk Mutect2 \
-R $fastaFile \
-I $inputFile \
--germline-resource $germlineResource \
-O $resultDir/$sampleName.Mutect.vcf.gz

Mutect2 uses a statistical model based on the available reads. Unlike with some other tools, there is not threshold based on depth or VAF or QUAL. The main statistic that does govern emission is the LOD (Log ODds) which is the log10 of the odds ("ratio") that the variation being inspected is a tumor relative to the probability that it is a technical artifact. The default threshold for this is 3.0 (corresponding to a odds ratio of 1:1000).
Perhaps you'd like to see lower probability results, in which case you should lower the threshold:
Using --tumor-lod-to-emit 2.0 should give you many more results, many of them might be wrong. The LOD isn't well calibrated, so I'd take it with a grain of salt....validate the results you get.


text:
How many false positive duplicates are marked using just the position of first unclipped base?
In the popular picard MarkDuplicates tool, a read is marked as a duplicate if it has the same position as another read starting from their first unclipped base in the 5'-direction. The implication is that if two reads have this property, it must be much more likely they arose from some "artificial" duplicates (PCR, optical, etc.) rather than due to sheer coincidence that two independent fragments of DNA happened to have this property and were sequenced in the same run.
Does anyone have links to studies or other data that attempts to measure how often this mechanism for labeling duplicates fails? I.e. how frequently does a read get marked a duplicate in this way "incorrectly," meaning it was actually just a coincidence and really did correspond to an independent fragment that should be considered in your analysis.
I imagine this answer can vary quite a bit based on the sequencing conditions: library prep, machine, etc. I was hoping for at least a ballpark estimate in the case of some "standard" modern short read sequencing methods. But it should be possible to study this in particular cases, e.g. using UMIs, doing PCR vs non-PCR runs, etc.

There's a tool in Picard (full disclaimer: I wrote it...) called CollectIndependentReplicateMetrics which performs the test you are looking for, I think.
It measures the percent of duplication that is "independent", i.e. which arose from randomly sampled original fragments that happen to be duplicates, rather than due to technical (optical, PCR, bridge-amp...) reasons.
The way it works is by looking at duplicate-sets of size 2 (and separately, 3) over heterozygous sites in the bam (which is assumed to be a diploid sample) it classifies the sets into "homogeneous" and "heterogeneous" based on the alleles of the hetrozygous variant. From the relative sizes of these two classes, one can work out the percent of duplication that is "independent" since independent duplication will give rise to heterogeneous sets (sometimes), while dependent duplication will never do that.
The assumptions are quite simple actually:

a technical duplicate will have the same allele
an independent replicate will have the same allele 50% of the time

So the rate of independence is 2x the rate of heterozygosity, within duplicate sets of size 2.
For sets of size 3 the math is a little more complex, so I wrote it up and docuemnted it.
Here's the verbatim explanation I wrote at the time (in the source-code of the Metric file...):
Explanation of calculation of independent replication rate from heterozygosity rate (within triplicate sets):

We assume the there are two types of replication events:

- those that are "independent", such that we just happen to get 2 fragments from the exact same region
(These get a random allele so effectively change the allele with probability 0.5), and
- those that are "artifactual" (do not change the allele)

We skip sets that have unexpected alleles as they do not fit our model. In the following we use the term "duplicates" to indicate that
the read-pairs would be marked as duplicates, not that they actually are technical duplicates.

To reach a triplicate set, 2 replication events are required so there are the following options, assuming
that the independent replication rate is x (thus the artifactual is 1-x). We assume a diploid organism with no bias towards either allele
in heterozygous sites, so an idependent replication will result in the other alleles in half the cases:

0 -> 0, 1 happens with probability x/2, therefore
0 -> 0, 0 happens with probability 1-x/2

(This is the explanation that is required for calculating the independent replication rate from doubleton sets...quite simpler)

Without loss of generality we assume that we "start" with allele 0 and that 1 is the other allele.

Each of the resulting alleles ("First" or "second") can replicate (each with probability 0.5) so we get:

from first row:
0=>1 1 (0 replicate to a 1) with probability  x/2 * 0.5 * x/2
0=>0 1 (0 replicate to a 0) with probability  x/2 * 0.5 * (1-x/2)
=====subtotal = x/4

0 1=>0 (1 replicate to a 0) with probability x/2 * 0.5 * x/2
0 1=>1 (1 replicate to a 1) with probability x/2 * 0.5 * (1-x/2)
=====subtotal = x/4

from second row:
0=>1 0 (first 0 replicate to a 1) with probability (1-x/2) * 0.5 * x/2
0=>0 0 (first 0 replicate to a 0) with probability (1-x/2) * 0.5 * (1-x/2)   <======= Homogeneous set
=====subtotal = (1-x/2)/2

0 0=>1 (second 0 replicate to a 1) with probability (1-x/2) * 0.5 * x/2
0 0=>0 (second 0 replicate to a 0) with probability (1-x/2) * 0.5 * (1-x/2)  <======= Homogeneous set
=====subtotal = (1-x/2)/2

total is x/2 (from first two sub-totals) + (1-x/2) (from last two sub-totals) = 1

We differentiate between a heterogeneous result (with 0's and 1's in the set) and homogeneous results
(all zeros, since we assumed WLOG that we start with 0)
The probability of a homogeneous set is therefore the sum of the two only homogeneous results

P(hom) = (1-x/2) * 0.5 * (1-x/2) + (1-x/2) * 0.5 * (1-x/2) = (1-x/2)^2 = y

where y = P(hom) is the rate of homogeneity within triplicate sets.

Solving for x we find that 2 * ( 1 - sqrt(y) ) = x.

I should also mention that this tool is based on the DuplicateSetIterator which is a really convenient tool (from htsjdk) for exploring duplication.

Hi Yossi, thanks for this! I just tried running this tool on a public sample using an Illumina 30x bam we have on a couple random chromosomes (14, 18, 20; 3 segfaulted for some reason). I ended up seeing biSiteHeterogeneityRate in the range 0.000544 to 0.000751 and triSiteHeterogeneityRate 0.001005 to 0.001262. If I'm understanding correctly, does this suggest (if my one sample is representative) that for this particular sequencing process we'd expect something like 0.05%-0.13% false positive duplication rate (so around 1/1000 duplicate reads)?

I'm assuming the tool has some built-in prior distribution that genuine duplicate-sets of size 2 at a het site should be 100% either ref or hom, but then tracks the rate at which the duplicates disagree at those sites? Does the tool have some other fancy techniques for trying to disambiguate the cases where you by chance have false positive duplicates but they both agree on the allele? I suppose one could take the "FP rate" from hets and extrapolate a similar rate to these sites and estimate some sort of overall FP duplication rate...

I documented the "model" in the tool itself...let me add to my answer.

@Ricky the numbers in the metric file that yoyu need are independentReplicationRateFromTriDups and independentReplicationRateFromBiDups. The tell you the estimated fraction of duplication which is not really a duplicate. I.e it came from a different original molecule.

As to whether one can recover the false positives, sure...but only over heterozygous sites, and only about half of them, as they would have to have a different allele from their sister reads. This also means that you already know what the genotypes are...so the value of extracting these reads is dubious....

The "heterogeneity" metric only tells you what fraction of the duplicate sets of size 2 & 3 were "heterozygous", i.e. didn't have reads that agreed on the allele. it's a step in the direction of getting the independent replication rate, but not the whole story (see explanation in the answer)

Thanks for clarifying @Yossi and modifying your answer. This is great. I'll have to think through the stats a bit more to get a final estimate on the FP rate from what I have (after I'm done traveling), but this tool is a nice convenience to have when thinking about this problem.

Sure thing @Ricky, and if you have any ideas of how to make the tool more accessible, or useful, I'm all ears....


text:
How to calculate pharmacophore features?
I asked this question already in the matter modeling forum (here). But i hope this is a better place because the topic is more biological.
I want to use the enzyme ß-Glucuronidase with the pharmit database. I uploaded my file for the receptor from pdb file. How can I get the pharmacophore features? So one thing is the enzyme is quiet popular, so is there any database with the features? I tried it with a publication but the given link is not working. It checked the site from the author but I did not found the database there.
And if there is no database or maybe there is one but the enzyme is not included, how can I get the pharmacophore features for my calculations?

Just put a know ligand inside the active site of your pdb , the webserver will calculate the Features of the ligand/pharmacophore then you can select the ones you want and the database where to search for molecules with the same Features. All the software used to crate a Features file seems to be not for free or really old (not sure if you can trust that code anymore) if you need some particular Features , you'll have to draw a FAKE LIGAND by yourself that carries them ... but thats its just my naive idea – –
se my resul using E.coli 3LPF pdb carrying Z77: 1-[(6,7-dimethyl-2-oxo-1,2-dihydroquinolin-3-yl)methyl]-1-(2-hydroxyethyl)-3-(3-methoxyphenyl)thiourea ligand and selecting first 4 Pharmacophire Features extracted by the webserver from it and running them against the CHEMBL32 database, first hit selected.

EDITED
The instructions above didnt work for me , loading 3lpf.pdb from my PC to the webserver for some reasons doesn't get the Z77 ligand. So the only option is to draw the Features by yourself using the Features menu, played a little bit with it but not well enough to get results , also the tricky part here would be to set the right coordinates for your ligand-Features, but overall the Tool looks very good, really fast and with a nice intuitive interface.
Still would like to get ligands from PC uploaded pdb file too , so that I could test my Fake-Ligand theory.

RE-EDITED :
found a third way :
load your pdb from your pc , save the session , modify the pharmit.json file you save adding the FAKE-LIGAND. To have an idea of how to put your FAKE-LIGAND inside the json file , have a look at  a saved session from a session started suppliying a pdb code instead of a file. In my experiment I loaded the modified Z77 ligand stiched to a session started from a downloaded 3lpf.cif and saved to a .json

auchhh !!!!! it does load the ligand but doesn't create the Features , don't know if it is because the ligand is cut in 2 pieces or because it just doesnit work like when I uploaded a .pdb containing a ligand.
At least you have a way to position the Features you'll add by hands afterwards.
***think that this one wont never works because the Features for the ligand are saved in the .json (first entry in the dictionary) the webserver calculates them once and then save them to be reloaded); yes it does do that , I can delete the ligand entry from .json and reload the session getting the Features but not the ligaand.
FOURTH EDIT (last one I swear) :
Tried new approach: loading Features menu uploading Z77 from https://www.rcsb.org/ligand/Z77 page, both model and ideal molecule as .sdf file. The ligand gets uploaded and Features calculated by the webserver , obviously the model molecule doesn't get located in the active site of the molecule, so I guess you could draw yor FAKE ligand using something like PyMOL exactly where you want it to be and then upload it with the load Features menu.
My uploaded Z77_ideal.sdf pic:

So that would provide ground for my FAKE-Ligand approach


text:
Variant allele frequency versus population frequency
In gnomAD related article
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9160216/#:~:text=Allele%20frequency%20is%20calculated%20by,can%20differ%20substantially%20between%20positions
I road this section

Allele frequency is calculated by dividing allele count by allele
number, hence, allele frequency represents the frequency of
confidently sequenced haplotypes that carry the allele in question
(because coverage and sequencing quality varies across the genome, the
allele number can differ substantially between positions). Allele
frequency is not equivalent to the percentage of individuals that
carry the allele, but is a suitable value for expressing the frequency
of a variant in the general population. The number of individuals
carrying a variant will depend on the number of heterozygous and
homozygous individuals but can be calculated from the data provided in
the variant table.

Does this mean, for getting the frequency of each allele in gnomAD, one could use Allele frequency (AF) column or must calculate population frequency?

What VCF fields are you comparing? I see the gnomAD has the overall AF field and various AF fields calculated by sample subset.

How did you subset gnomAD VCF to non-finnish samples only?

You might find gnomAD's non-cancer AF fields useful. I'll try to decipher the formula. AN/2 is the number of samples (assuming all samples are diploid). AC - num_homo gives you 2*(num(hets)). Dividing this by number of samples gives you frequency of samples with at least one ALT allele.

gnomAD provides a global Allele Frequency as well as AF values for most subsets. gnomAD does not provide individual genotypes (as far as I know), so calculating custom AF for non-precomputed subsets might not be possible.
Since you're interested in the NFE (non-Finnish European) sub-population, the nfe_*_AF fields might be of interest to you.

As far as calculating the fraction of non-WT samples from num_homozygous_samples, AC, and AN goes, your formula seems spot on. Here's my deconstruction of the formula:
Your formula is: (AC - num_homozygous_samples)/(AN / 2))
AN/2 is the number of samples (assuming all samples are diploid). AC - num_homo gives you 2 * num_hets. Dividing this by number of samples gives you frequency of samples with at least one ALT allele.
For example, let's take a dataset of 200 diploid samples where 100 are WT, 75 are het and 25 are hom-alt for a biallelic variant.
Since there are 200 diploid samples, AN is 400.
AC for the ALT allele is (75 * 1) + (25 * 2) = 125.
Since 100 are WT, it is known that fraction of samples with mutation = 1 - (100/200) = 0.5. This number can be arrived using your formula:
(AC - num_homozygous_samples) / (AN/2)
 = ((125 - 25) / (400/2))
 = 100 / 200
 = 0.5


text:
Error in GSEA expresion data set file
I'm want to perform an GSEA of some RNA-seq data but when I upload the expression data set file this problem appears but I don't know how to solve it.
Thanks for your help.
Pol
Error Details
---- Full Error Message ----
There were errors: ERROR(S) #:1
Parsing trouble
java.lang.NumberFormatException: ...

---- Stack Trace ----
# of exceptions: 1
------For input string: "description"------
java.lang.NumberFormatException: For input string: "description"
    at java.base/jdk.internal.math.FloatingDecimal.readJavaFormatString(Unknown Source)
    at java.base/jdk.internal.math.FloatingDecimal.parseFloat(Unknown Source)
    at java.base/java.lang.Float.parseFloat(Unknown Source)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.AbstractParser.parseStringToFloat(AbstractParser.java:250)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.AbstractParser.parseFieldsIntoFloatArray(AbstractParser.java:361)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.AbstractParser.parseTextMatrixToDataset(AbstractParser.java:277)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.TxtDatasetParser.parse(TxtDatasetParser.java:108)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.readDatasetTXT(ParserFactory.java:200)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:758)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:735)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserWorker.doInBackground(ParserWorker.java:53)
    at java.desktop/javax.swing.SwingWorker$1.call(Unknown Source)
at java.base/java.util.concurrent.FutureTask.run(Unknown Source)
at java.desktop/javax.swing.SwingWorker.run(Unknown Source)
at java.base/java.util.concurrent.ThreadPoolExecutor.runWorker(Unknown Source)
at java.base/java.util.concurrent.ThreadPoolExecutor$Worker.run(Unknown Source)
    at java.base/java.lang.Thread.run(Unknown Source)

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.

Please show the first few lines of your input files. Given this is a file format error problem, it will help to know the current format of your files.


text:
What does a question mark ("?") mean in Picard metrics files when I expect a number (integer, float, etc.)
Does anyone know what conditions must be true for Picard to put a ? where I expect a number?
Here is an example from the output of the tool CollectWgsMetrics:
❯ cat sample.wgs_metrics.txt | head -n8 | tail -n2 | verticalize | grep FOLD
FOLD_80_BASE_PENALTY    1.33476
FOLD_90_BASE_PENALTY    2.288159
FOLD_95_BASE_PENALTY    ?

Code searching (link) didn't reveal anything, and googling isn't returning anything obvious either.

[Verticalize is a little CLI tool I've forever been using by Pierre Lindenbaum: verticalize. Easy to compile and I toss it in my ~/bin every time I get a new computer.]

Picard (actually HtsJdk) uses a formatting class (via the MetricsFile class) that emits "?" for NaNs (and Infs).
In your case, there seems to be a division-by-zero error during the calculation of the metric (likely due to zero coverage at some small percentage of bases in your input). This resulting Nan is emitted as a "?".


text:
Converting .tsv to .meg for visualization
This is an add-on to the previous question I asked here.
I was able to deploy a "chunking" strategy to my kmer matrix and was able to get a file with the distance measurements! From what I could tell by using the head command was that the data looks like what it did before when I was using this approach on a smaller dataset so now I'm going to use this data to generate a tree using MEGAX.
Normally when I'd convert the tsv file to a meg file I'd manually add in the meg headers/requirements for conversion. However, I can't open the distance tsv file (its ~23 GB) to make these changes.
Is there an easy solution to this problem? How do you convert a tsv file to a meg file without opening the file?
EDIT:
For more context here is the format for each file.
Currently, my tsv file looks like this:
0 1 2 3 4 5 6 7 8 9 10 ... 178
1 0 # # # # # # # # #  ... #
2 # 0 # # # # # # # #  ... #
3 # # 0 # # # # # # #  ... #
4 # # # 0 # # # # # #  ... #
5 # # # # 0 # # # # #  ... #
6 # # # # # 0 # # # #  ... #
7 # # # # # # 0 # # #  ... #
8 # # # # # # # 0 # #  ... #
9 # # # # # # # # 0 #  ... #
10 # # # # # # # # # 0 ... #
...
178 # # # # # # # # # #  ... 0

This is a jaccard distance matrix of 178 samples. The # represent the distance values for each sample compared to the other. I was able to calcualte the distance using pandas and scipy: matrix = pd.DataFrame(squareform(pdist(data, 'jaccard')))
And I would like it to look like this:
#mega
TITLE: Jaccard distance of 178 samples

sample_1
sample_2
sample_3
...
sample_178

distance_1
distance_1 distance_2
distance_1 distance_2 distance_3
distance_1 distance_2 distance_3 distance_4
...
distance_1 distance_2 distance_3 distance_4...distance_178

I added some extra information to give more context to this issue! Let me know if you'd like me to clarify something further or if you'd like more information. I appreciate the help!

Thanks @rimo .. okay this ain't the data set I imagined. I'll need to reproduce it. Are you using SciPy, i.e. where is the pdist() method coming from?

Note, its a while since I've used scipy and I don't fancy loading an deleting for this test ... so in the absence of pdist() the solution would be ...
import numpy as np
import pandas as pd
from pathlib import Path

data = [[1,2,3],[0,5,6],[7,0,9],[10,11,0]]
arr = np.tril([data]).astype('float')
arr[arr == 0] = 'nan'
arr = np.where(np.isnan(arr), '', arr)
columns = ['a','b','c', 'd']
file = Path('/Volumes/Data2/test/test1.meg')
if file.is_file():
    file.unlink()
with open(file, 'a') as fout:
    fout.write("#MEGA\n")
    for col in columns:
        fout.write(col + "\n")
    fout.write("\n0\n")
    np.savetxt(fout, arr.reshape(3,-1), fmt='%s')

output
#MEGA
a
b
c
d

0
1.0 0.0
5.0 7.0 0.0
9.0 10.0 11.0 0.0

That works. Keep in mind where it says arr.reshape(3,-1) the 3 should be the number of columns -1 in your matrix: in your case this appears to be 177 (178 - 1).
The gist is to use np.tril, which doesn't exists in pandas. Hence don't use pandas.
At a guess with scipy then it would be ...
import numpy as np
import pandas as pd
from pathlib import Path
from scipy import pdist

tmp = [[1,2,3],[4,5,6],[7,8,9],[10,11,12]]

data = pdist(tmp, 'jaccard')
arr = np.tril([data]).astype('float')
arr[arr == 0] = 'nan'
arr = np.where(np.isnan(arr), '', arr)
columns = ['a','b','c', 'd']
file =Path('/Volumes/Data2/test/test1.meg')
if file.is_file():
    file.unlink()
with open(file, 'a') as fout:
    fout.write("#MEGA\n")
    for col in columns:
        fout.write(col + "\n")
    fout.write("\n0\n")
    np.savetxt(fout, arr.reshape(3,-1), fmt='%s')

You might need a quick perl -p -i -e 's/[" "]+/ /g' test1.meg or even perl -p -i -e 's/[" "]+/\t/g' test1.meg to set the file for MEGA. I dunno why but \s or " " in the transform part of perl pie isn't working for me.


text:
Chromosome position -> ENSP
I'd like a tabix file that, for every position in the genome, gives me a list of every ENSP within 5kb of that position.
This should be easy right... I know VEP /kinda/ gives me this, but I'd rather not run it for every postition in the genome
Can I use biomart to get the start/end position for every ENSP, then bedtools to somehow expand those coords?

This seems like an XY problem. What is the overarching reason for this task?

What's an "ENSP"? Do you mean Ensembl protein IDs? What species are you working with? Which genome version? What you describe is simply every single annotated protein, so why aren't you just taking all proteins? Please [edit] your question and explain what the end goal is so we can understand what you need better.

I'm assuming by "ENSP" you mean Ensembl proteing (which have codes starting with ENSP, then a number). Is that correct?

I used BioMart (and then tabix) to get the answer to this question:
if (!requireNamespace("BiocManager", quietly = TRUE)) {
  install.packages("BiocManager")
}

BiocManager::install("biomaRt")

## F**kery:
install.packages("devtools")
devtools::install_version("dbplyr", version = "2.3.4")

library(biomaRt)

listEnsembl()

ensembl <- useEnsembl(
  biomart = "ensembl",
  dataset = "hsapiens_gene_ensembl",
  GRCh = 37
)

chrs <- c(1:22, "X", "Y", "MT")

genes <-
  getBM(
    attributes = c(
      "chromosome_name",
      "start_position",
      "end_position",
      "ensembl_gene_id",
      "ensembl_transcript_id",
      "ensembl_peptide_id",
      "hgnc_symbol"
    ),
    filters = "chromosome_name",
    values = chrs,
    mart = ensembl
  )


text:
Calculating the segregating sites via tskit versus counting
Can you help me understand why I am getting different values when I calculate the segregating sites from tskit api versus manually counting? Thanks!
Example data generation:
import sys, tskit, msprime

ts = msprime.sim_ancestry(samples=3, ploidy=2, sequence_length=10, random_seed=2)
ts = msprime.sim_mutations(ts, rate=0.1, random_seed=2)

Example results: I get 5 from API and 6 from manual counting.
>>> ts.segregating_sites(span_normalise=False)
5.0
>>> n=0
>>> for tree in ts.trees():
...     for site in tree.sites():
...         for mutation in site.mutations:
...             node = ts.node(mutation.node)
...             pop_id = node.population
...             n+=1
...
>>> n
6


text:
Q8 Accuracy Evaluation metric mathematical expression used in Protein Secondary Structure Prediction
I am working on Protein Secondary Structure Prediction using Recurrent Neural Network (GRU) model. I came across few open source projects which already worked on the similar problem. All of them are using Q8 Accuracy as evaluation metrics and the custom code they used for calculating Q8 accuracy is given below.
def accuracy(y_true, y_pred):
  y = tf.argmax(y_true, axis =- 1)
  y_ = tf.argmax(y_pred, axis =- 1)
  mask = tf.greater(y, 0)
  return K.cast(K.equal(tf.boolean_mask(y, mask), tf.boolean_mask(y_, mask)), K.floatx())

I would like to know the mathematical expression of Q8 accuracy corresponding to the above code snippet. The above code source project reference https://github.com/idrori/cu-ssp


text:
Discrepancy between UniProt and AlphaFold data for MlaD domains – Seeking guidance on automating entry sorting based on accurate domain count
Hey fellow researchers,
I recently conducted an analysis on the UniProt database, specifically focusing on entries with MlaD domains. After downloading the data and opening it in an Excel file, I noticed a significant discrepancy with the AlphaFold predictions.
For instance, entry "A0A336N986" in UniProt indicates only 1 MlaD domain (link: UniProt Entry), but when cross-referencing with AlphaFold (link: AlphaFold Entry), it reveals 5 MlaD domains with high confidence values.
Considering the large dataset with over 100k entries, manually sorting them based on the accurate count of MlaD domains seems impractical. I'm reaching out to the community for guidance on automating this process. Has anyone encountered a similar issue, or does anyone have suggestions on tools, scripts, or strategies to efficiently sort UniProt entries based on the correct number of MlaD domains?
Any insights or advice would be greatly appreciated! Thanks in advance.


text:
Population Allele Frequency
In UK biobank research environment
In GENOMICS tab of UK BRAP, if I type 7:140476600-140480852 coordinate in search bar and the same time after filtering my cohort for this coordinate, in ALLELE TABLE I see columns like:
Cohort AF Cohort Allele Count Population AF Population Allele Count GnomAD AF

I Googled a lot but no clue where population allele count and frequency come from
Does anyone have any experience with this?

I don't know specifically about the UKBB RAP, however, Population Allele Count is the count of that specific allele genotyped in the population in question. In this case, it most likely refers to the population that your cohort is part of. The population allele frequency is the same thing but noted in terms of frequency, i.e. divided by 2 times the number of genotyped individuals in that population.


text:
Dealing with X Residues in FASTA files prior to folding analysis
We are working with a large dataset of proteins for folding/docking tests across various tools (AlphaFold, ESM2, RoseTTAFold). In some of the FASTA files, there is an X for a non-standard residue. An example of this is in the 6SEN_1 FASTA file below. Is there a standard method of dealing with these non-standard residues prior to analysis to minimize impact on the folding/docking results?
>6SEN_1|Chains A, B|Transcriptional enhancer factor TEF-3|Homo sapiens (9606)
GPRSVASSKLWMLEFSAFLEQQQDPDTYNKHLFVHIGQSSPSYSDPYLEAVDIRQIYDKFPEKKGGLKDLFERGPSNAFFLVKFWADLNTNIEDEGSSFYGVSSQYESPENMIITCSTKVCSFGKQVVE**X**VETEYARYENGHYSYRIHRSPLCEYMINFIHKLKHLPEKYMMNSVLENFTILQVVTNRDTQETLLCIAYVFEVSASEHGAQHHIYRLVKE

X  Any - Uncommon or Unspecified

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8942982/   ---> The role of lysine palmitoylation/myristoylation in the function of the TEAD transcription factors , Just change it to a LYS ??'

Similar to @pippo1980 comment: start with a related standard amino acid, fold it however you are doing that, and then mutate to the nonstandard in the folded structure. You can do this with a toolkit like pdbfixer for standard AAs, you might need to do something like this to directly add the atoms for nonstandards, using an example pdb structure of the nonstandard AA.
You would then need to follow that with a round of energy minimization in the structure to resolve any large problems in the structure, using a toolkit like OpenMM.
If X truly means that you have no idea what's there, then... that's harder. In this specific case lysine may be a best bet, but it may be hard to generalize between e.g. methyllysine, phosphotyrosine, phosphoserine, etc. They're pretty different structurally, so having a good idea of what family of residue you would expect would be helpful.
The safest thing to do would be to go back to the transcript and in vitro translate to get an idea of what that X is. Though that's getting a little complicated compared to just replacing it with a leucine or whatever.


text:
Are there any simple ways to call copy number variants from a Nebula Genomics WGS data?
I want to find my copy number variants. I got DNA testing from the Nebula Genomics and received the raw datas in .cram, .fastq, .vcf formats. Are there any simple and easy method to call my copy number variants from one of these files? I tried GATK, but I was unable to collect Read Counts from 1000 genomes datas.

Yes and no. Yes, there are tools that can do this, but no, they do not work very well. Short read NGS sequencing just isn't the right tool to find large structural variants like CNVs. As a rule of thumb, you can expect around 30% accuracy when compared to other experimental approaches such as MLPA or array CHG. That said, yes, there are indeed various tools that can do this. My personal favorite is delly because it is accurate, actively maintained and since it is written in C it is quite fast.
Once you have downloaded and installed it, you'll likely need to convert your cram to bam as I don't think it can deal with cram:
samtools view -b foo.cram > foo.bam

Then you can run delly on it with (assuming your cram was aligned against hg38):
delly call -x path/to/delly/excludeTemplates/human.hg38.excl.tsv \
   -g path/to/hg38.fa -o out.bcf foo.bam

You will also need to download a copy of the human genome fasta file that Nebula used to generate your cram. You can ask them where to get one, and change path/to/hg38.fa with the path to that file.


text:
Dog (Canis lupus familiaris) MT gene annotations
I have an RNA-Seq and scRNA-Seq dataset from dog samples and pseudo-aligned them to the dog reference transcriptome (Ensembl release 109). When filtering for mitochondrial gene percentage in scRNA-Seq I couldn't find the "MT" annotation in this new release. Does anyone have any information about this? Did they not include it or am I doing something wrong somewhere (the MT annotations are present in the CanFam3.1 genome built but not in the newer ones).
I could not find anything on the Ensembl release details or blog.
Thank you :)
This is a part of the code that I used to search for the mitochondrial genes:
searchFilterOptions(mart=ensembl_canis, filter="chromosome_name", pattern="MT")

Currently, there is no available mitochondrial sequence for the ROS_Cfam_1.0 assembly:
https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_014441545.1/

Thank you for your reply Ben!


text:
RNAseq and TTseq results
I'm working with two sets of data. One of two sets of data, one contains the results of a RNAseq analysis and the other contains the result of TTseq analysis. I want to check which genes overlap in both data sets to check if the are affected at stabilization level or at a transcriptional level. If anyone knows how the coding pipeline it would help me. I work on linux.

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:
OsError in virtual box snakemake pipeline for array-analysis
The rule is written for calling genotypes using array analysis cli for illumina gsa chips.
import configparser
config=configparser.ConfigParser()
config.read(".//config")
homepath=config.get('Paths','home_dir')
inputpath=config.get('Paths','input_dir')
outputpath=config.get('Paths','output_dir')
Sentrix_ids=config.items('Sentrix_ids')
sentrixids=[]
positions=['R01C01','R01C02',
                        'R02C01','R02C02',
                        'R03C01','R03C02',
                        'R04C01','R04C02',
                        'R05C01','R05C02',
                        'R06C01','R06C02',
                        'R07C01','R07C02',
                        'R08C01','R08C02',
                        'R09C01','R09C02',
                        'R10C01','R10C02',
                        'R11C01','R11C02',
                        'R12C01','R12C02']

rule all:
        input:
                expand(outputpath+"/{sids}_{pos}.gtc",sids=sentrixids,pos=positions)

rule idat2vcf:
        input:
                idatpath=homepath+"/ss.csv",
                clusterfile=homepath+"/input/"+"GSA-24v3-0_A1_ClusterFile.egt",
          manifestfile=homepath+"/input/"+"GSA-24v3-0_A1.bpm"

        output:
                outputfol=outputpath+"/"+"gtcfiles"
        shell:
                "bash /scripts/idattogtc.sh {input.idatpath} {input.clusterfile} {input.manifestfile} {output.outputfol}"

#!/usr/bin/bash
exec_path="/home/user14/Downloads/array-analysis-cli-linux-x64-v2.1.0/array-
analysis-cli/array-analysis-cli"
idat_path=$1
cluster_file=$2
manifest_file=$3
output_fol=$4
${exec_path} genotype call \
  --bpm-manifest ${manifest_file} \
  --cluster-file ${cluster_file} \
  --idat-sample-sheet ${idat_path} \
  --output-folder ${output_fol}

When I run snakemake, it gives the following error:
Traceback (most recent call last):
  File "/home/user14/miniconda3/envs/snakemake/lib/python3.6/site-packages/snakemake/__init__.py", line 469, in snakemake
    force_use_threads=use_threads)
  File "/home/user14/miniconda3/envs/snakemake/lib/python3.6/site-packages/snakemake/workflow.py", line 447, in execute
    self.persistence.cleanup_shadow()
  File "/home/user14/miniconda3/envs/snakemake/lib/python3.6/site-packages/snakemake/persistence.py", line 128, in cleanup_shadow
    shutil.rmtree(self.shadow_path)
  File "/home/user14/miniconda3/envs/snakemake/lib/python3.6/shutil.py", line 490, in rmtree
    onerror(os.rmdir, path, sys.exc_info())
  File "/home/user14/miniconda3/envs/snakemake/lib/python3.6/shutil.py", line 488, in rmtree
    os.rmdir(path)
OSError: [Errno 26] Text file busy: '/media/sf_vm_centfiles/workflow/.snakemake/shadow'

Are you running Snakemake within a shared folder in Virtual Box by any chance?  Here's somebody with what sounds like a similar problem.


text:
Does the `mmcif/pdb` format store atomic coordinates in units of Ångstrom, nanometers or something else?
I was surprised as to how hard it was to find a reference on this among the docs.
Does the Protein Data Bank mmcif/pdb-x archival format use Angstrom to specify cartesian x,y,z coordinates of atoms or something else or abstract units (perhaps defined per each mmcif _data frame?
How is the origin of the grid calculated -- for crystallographers out there?
Thank you

From google I got this as first hit :
PDB_format_1992.pdf :

Hopefully they keept the same Unit of measurement during the data/file format evolution/extension
It seems the case for mmcif Data Item _atom_site.Cartn_x :

Some history here :
The Protein Data Bank: a historical perspective , The Protein Data Bank
AND :
PY  - 1971
DA  - 1971/10/01
TI  - Crystallography: Protein Data Bank
JO  - Nature New Biology
SP  - 223
EP  - 223
VL  - 233
IS  - 42
SN  - 2058-1092
UR  - https://doi.org/10.1038/newbio233223b0
DO  - 10.1038/newbio233223b0
ID  - ref1
ER  -

About the axis origin CHARMM-GUI has a nice introductory lesson to PDB format (see CRYTST1 Record Section) Lesson 1: Introduction to PDB Format
and stolen from the web, free part of an old book :

Thanks, pip. Missed the "angstroms" part when browsing through the spec. The paper and book references are very much appreciated.


text:
Recommendations on Motif scoring functions
I am searching for specific transcription factor binding locations on DNA, and ranking them according the their scores.
For this, I am in search of a tool or method that can generate motif scores for arbitrary-length DNA sequences (denoted as 'x') using a scoring function Motif_score(x, Motif_M_PWM). The purpose is to calculate motif scores in a moving window fashion along the entire DNA sequence. The 'Motif_M_PWM' parameter typically represents the probability weight matrix (PWM) associated with a specific motif 'M.'
One common method to calculate the motif score for a given short sequence 'x' (which has the same length, n, as the motif itself) is by using the formula:
$motif\_score = \sum_{k=1}^{n} \log\frac{PWM\_Probability(x(k), k)} {0.25}$
In this formula, 'x(k)' represents a specific nucleotide in sequence 'x,' and 'k' denotes the position of the nucleotide in the  PWM and x. The '0.25' term accounts for the random chance of any one of the four nucleotides (A, C, G, T) occurring. I think this scoring method is commonly used by tools like HOMER.
I am curious to know if there are alternative scoring functions or methods (tools besides HOMER) available for calculating motif scores on DNA sequences.


text:
Concatenate multiple fasta files into individual files per unique entry
I have several splitted fasta files like this:
TK273759.fastaaa
TK273759.fastaab
TK273759.fastaac
TK273759.fastaad
TK273759.fastaae

I want to concatenate these fasta files into one file like this:
TK273759.fasta

What command should I use to make this happen? I appreciate your help.

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.

Please try searching online before asking for help. A simple cat will do the trick here:
cat TK273759.fasta[a-z][a-z] > TK273759.fasta

Thank you for your help!

Please accept the answer to provide closure to the post.


text:
Defining out-groups in MEGA
I am using the MEGA software for trimming protein sequences and building phylogenetic trees.
How can I define out-groups of phylogenetic tree in MEGA? Is there a way to calculate the distances of each species from the root?

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.

Hmmm.... the issue is that the precise location of the root - unless it's a molecular clock - is arbitrarily assigned within limits. You can designate populations ('clades') and ask for the genetic distance between the clades, it's one of the menus header options. To do what you want you need to build the tree (neighbor-joining method) and assess it in the pop-up window. Switch the view of the tree to 'genetic distance' and then manually add the branch lengths for each root - tip distance. One of the Beast supporting packages use to calculate the root tip distance, I think this option was removed in the upgrades.