text:
Gencode PolyA feature annotation GTF gene_id not the same as Comprehensive gene annotation GTF
I have downloaded the Gencode PolyA feature annotation GTF here. The first 10 lines are like this:
##description: evidence-based annotation of the human genome (GRCh38), version 41 (Ensembl 107) - polyA features
##provider: GENCODE
##contact: gencode-help@ebi.ac.uk
##format: gtf
##date: 2022-05-12
chr1    HAVANA  pseudo_polyA    134957  134962  .   -   .   gene_id "383460"; transcript_id "383460"; gene_type "pseudo_polya"; gene_name "383460"; transcript_type "pseudo_polya"; transcript_name "383460"; level 2;
chr1    HAVANA  polyA_site  141474  141475  .   -   .   gene_id "383461"; transcript_id "383461"; gene_type "polya_site"; gene_name "383461"; transcript_type "polya_site"; transcript_name "383461"; level 2;
chr1    HAVANA  polyA_signal    141477  141482  .   -   .   gene_id "383462"; transcript_id "383462"; gene_type "polya_signal"; gene_name "383462"; transcript_type "polya_signal"; transcript_name "383462"; level 2;
chr1    HAVANA  pseudo_polyA    184934  184939  .   -   .   gene_id "383463"; transcript_id "383463"; gene_type "pseudo_polya"; gene_name "383463"; transcript_type "pseudo_polya"; transcript_name "383463"; level 2;
chr1    HAVANA  pseudo_polyA    258530  258535  .   -   .   gene_id "383464"; transcript_id "383464"; gene_type "pseudo_polya"; gene_name "383464"; transcript_type "pseudo_polya"; transcript_name "383464"; level 2;

In this GTF file, what does the gene_id mean? In the whole genome GTF (same build) the gene_id is clearly the ENSEMBL gene_id:
##description: evidence-based annotation of the human genome (GRCh38), version 41 (Ensembl 107)
##provider: GENCODE
##contact: gencode-help@ebi.ac.uk
##format: gtf
##date: 2022-05-12
chr1    HAVANA  gene    11869   14409   .   +   .   gene_id "ENSG00000223972.5"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; level 2; hgnc_id "HGNC:37102"; havana_gene "OTTHUMG00000000961.2";
chr1    HAVANA  transcript  11869   14409   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_name "DDX11L1-202"; level 2; transcript_support_level "1"; hgnc_id "HGNC:37102"; tag "basic"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
chr1    HAVANA  exon    11869   12227   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_name "DDX11L1-202"; exon_number 1; exon_id "ENSE00002234944.1"; level 2; transcript_support_level "1"; hgnc_id "HGNC:37102"; tag "basic"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
chr1    HAVANA  exon    12613   12721   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_name "DDX11L1-202"; exon_number 2; exon_id "ENSE00003582793.1"; level 2; transcript_support_level "1"; hgnc_id "HGNC:37102"; tag "basic"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
chr1    HAVANA  exon    13221   14409   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_name "DDX11L1-202"; exon_number 3; exon_id "ENSE00002312635.1"; level 2; transcript_support_level "1"; hgnc_id "HGNC:37102"; tag "basic"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";

What I want to do eventually is to assign the PolyA sites to ENSEMBL gene_ids.

According to the GENCODE ftp README file, the value of the 'gene_id', 'transcript_id', 'gene_name' and 'transcript name' fields corresponds to a random identifier. The polyA features are not directly associated to any gene or transcript when manually annotated by Havana.
There is other file in the GENCODE release that includes links between polyA features and transcripts:
https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_41/gencode.v41.metadata.PolyA_feature.gz

ENST00000485748.5       3232    3233    chr1    1496201 1496202 +       polya_site
ENST00000485748.5       3211    3216    chr1    1496180 1496185 +       polya_signal
ENST00000474481.1       3367    3368    chr1    1496201 1496202 +       polya_site
ENST00000474481.1       3346    3351    chr1    1496180 1496185 +       polya_signal
ENST00000308647.8       2114    2115    chr1    1496201 1496202 +       polya_site
ENST00000308647.8       2093    2098    chr1    1496180 1496185 +       polya_signal

22. gencode.vX.metadata.PolyA_feature.gz:
  Manually annotated polyA feature overlapping the transcript 3'-end.
  1 - transcript id
  2 - transcript-based start coordinate of the polyA feature
  3 - transcript-based end coordinate of the polyA feature
  4 - polyA feature chromosome
  5 - polyA feature start coordinate
  6 - polyA feature end coordinate
  7 - polyA feature strand
  8 - polyA feature type ("polyA_site", "polyA_signal", "pseudo_polyA")


text:
Write a bash script to run gatk, fix errors with input, and rerun until completion
I have a bam file that I want to run through GATK's SplitNCigarReads tool. Because of the way the bam file was generated, the program will often fail, with an error message stating:

##### ERROR MESSAGE: Bad input: Cannot split this read (might be an empty section between Ns, for example 1N1D1N): ABCDEFG

My question is: how do I write a bash script to:

Start GATK on the input bam file.

When (and if) the error message arises, take the bad read (e.g., ABCDEFG) and use samtools to remove it,
samtools view -h input.bam | grep -v ABCDEFG | samtools view -bS -o intermediate1.bam

Index the new bam file

Run the same GATK command as #1, but with the newly generated bam file

Keep performing steps 1-4 until the GATK tool runs to completion

I realize this may be quite complicated and/or involved, so I don't expect someone to write the entire script for me. Are there any online resources dedicated to writing a script like this?
P.S. This is my first StackExchange question, please let me know if I did anything incorrectly or inefficiently.

Just wondering if you have an example BAM that produces the error that you could share with us? Alternatively, could you describe how you are generating your BAM files so that we can produce a file that will produce one or more errors when run?

Removing offending reads one by one and rerunning the tool every time will be very time-consuming, and strikes me as fundamentally the wrong approach: due to how BAM files are structured, removing individual reads is incredibly inefficient. Your time would be better spent devising a way to detect and remove the offending reads all at once.

may refer to page 217:  https://books.google.com.hk/books?id=JUp3EAAAQBAJ&pg=PA214&lpg=PA214&dq=bash+script+%22splitncigarreads%22&source=bl&ots=qf2xcKQL7Q&sig=ACfU3U1EQhIsOiGEnBndigXHEjnBkFX4aQ&hl=en&sa=X&ved=2ahUKEwi5t6XJyp36AhWQZd4KHaDUBDEQ6AF6BAgfEAM#v=onepage&q=bash%20script%20%22splitncigarreads%22&f=false

source: Stack Exchange

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

sort: votes


text:
Building a group representation of genes' ontologies
I have a list of genes (about 2000) from lines of drosophila subobscura, in the following format:

LOC117900589
LOC11788959
LOC11790331

I'm trying to represent this set of genes in a useful way.
So, I was recommended to represent the ontology function of these genes in a graphical way.
I was recommended this tool:
http://cbl-gorilla.cs.technion.ac.il/
One inputs a list of genes and it outputs this kind of graph:

The problem is that, this tool, GOrilla only accepts genes from drosophila melanogaster.
I have to either:

a) convert the genes id to a homologous gene in melanogaster
b) find an alternative tool

A solution to either a) or b) would be helpful. Or an alternative solution.
For a) I tried:

Ncbi's homologene: which doesn't seem to be working for a the genes ID I tried so far (I didn't try them all)
Blasting the gene sequence on melanogaster's genome and looking through the list manually to find an homologous gene.

So far these methods have been impractical for the number of genes I have.

Have you tried DAVID for functional annotation (Gene Ontology and KEGG terms)?
I usually use DAVID to extract GO and KEGG terms and plot the results using R.
If you also want to plot a protein-protein interaction (PPI) network for your gene list you can use STRING database to extract interactions and network and use Cytoscape for visualizing this result.

I tried putting a bunch of entrez id genes and I got nothing useful it seems? https://imgur.com/tZc6Vk8.png

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


text:
Aligning PacBio HiFi reads to reference genome using pbmm2
I am trying to align a yeast strain sequenced by PacBio HiFi reads to the reference genome S288C (https://www.ncbi.nlm.nih.gov/data-hub/genome/GCF_000146045.2/) using pbmm2 but for some reason I am getting an error that I don't understand.
I first indexed the reference sequence with the following code:
pbmm2 index S288Cseq.fa S288Cindx.mmi

And then I attempted to align using this line of code:
pbmm2 align S288Cseq.mmi 088.fastq.gz mapped088.aligned.bam

But I keep getting this error:
ERROR: unsupported extension on input file: S288Cindx.mmi
I get that error no matter if I use the reference file in either mmi or the fasta format. Can someone please help and let me know what I'm doing wrong? (I downloaded the tool using bioconda: conda install -c bioconda pbmm2 just in case that helps)

What happens if you try directly using the fasta in the alignment command? Indexing is optional.

Can you clarify what version of pbmm2 you are using? pbmm2 --version or something like that.

The same error comes up if I use the fasta file too. The version is:
pbmm2 0.8.1 (commit cbf547249)

I used minimap2 in its place for now and it worked that way but I want to see if I can still use this tool

Looking more closely, the index filename changes ( S288Cseq.mmi1 vs. S288Cindx.mmi). Maybe this is just a typo but it's worth confirming that nothing else is weird.


text:
Coordinate numbering between forward (+) and reverse (-) strand in Tombo
My question is related to the output of Tombo re-squiggling algorithm on nanopore data.
Given the position/coordinate on one strand (+ or -), I want to find the coordinate of its exact location on the opposite strand.
For example, let's say I want to find data/signal on specific coordinates on both strands on Tombo path after re-squiggling. Let's say position 13936 - 14036 on genome on + and - negative strands. I wrote the following code but the sequence obtained from Tombo is not similar at all. I want to know how I can find the coordinate corresponding to a given position on the opposite strand.
# one read from the + strand
mapped_start, mapped_end, read_id, fast5_path, _, _ = reads_per_win_pos.loc[rd_id, :]
>>> start
13936
>>> end
13942
>>> mapped_start, mapped_end, read_id, fast5_path, _, _ = reads_per_win_pos.loc[rd_id, :]
>>> mapped_start
11205
>>> mapped_end
20084
>>> f = h5py.File(fast5_path, 'r')
>>> tombo_events = list(f['Analyses/RawGenomeCorrected_000/BaseCalled_template/Events'])
>>> forward_seq = ''.join([e[4].decode() for e in tombo_events])
>>> start_idx = start - mapped_start
>>> start_idx
2731
>>> win_seq = seq[start_idx: start_idx + 100]
>>> forward_win_seq = forward_seq[start_idx: start_idx + 100]
>>> # one other overlapping read on the other strand (-)
...
>>> mapped_start, mapped_end, read_id, fast5_path, _, _ = reads_per_win_neg.loc[rd_id, :]
>>> mapped_start
10979
>>> mapped_end
20091
>>> f = h5py.File(fast5_path, 'r')
>>> tombo_events = list(f['Analyses/RawGenomeCorrected_000/BaseCalled_template/Events'])
>>> reverse_seq = ''.join([e[4].decode() for e in tombo_events])
>>> start_idx = start - mapped_start
>>> start_idx
2957
>>> reverse_win_seq = reverse_seq[start_idx: start_idx + 100]

Now the forward_win_seq and reverse (or even the reverse complement) of it are not similar:
>>> forward_win_seq
'GAGACGGAGCAGACCCATCTGCTACTGCCCTTTCTATAATAACTAAAGTTAGCTGCCCTGGACTATTCACCCCCTAGTCTCAATTTAAAAAGATCCCCAT'
>>> reverse_win_seq
'GGGCCAGGTGCCTGAGATCCATGTTCCATCCTACCTGCCTGACCTGCCCGGCATTGCCAACGACCTCATGTACAGTGCCGACCTGGGCCCCGGCATTGCC'
>>> complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
>>> com_seq = ''.join([complement[base] if base in complement.keys() else base for base in reverse_win_seq])
>>> com_seq[::-1]
'GGCAATGCCGGGGCCCAGGTCGGCACTGTACATGAGGTCGTTGGCAATGCCGGGCAGGTCAGGCAGGTAGGATGGAACATGGATCTCAGGCACCTGGCCC'

When Tombo re-squiggle and aligns a read it adds two attributes to it: mapped_start and mapped_end. Tombo also gives chromosomes and strands.
Imagine you are looking for some value of a signal[s] in position x to y on a chromosome on both strands. In this range, you probably have reads that span the positive strand (forward) and the reads that cover the negative strand (reverse).
Now let a given read that has mapped to a positive strand and covers the entire [x, y], have attributes: pos_mapped_start and pos_mapped_end. Similarly let the reads mapped to negative strands covering [x, y] have attributes: neg_mapped_start and neg_mapped_end.
For reads mapped to + strand:
start_idx = x- pos_mapped_start
end_idx = y - pos_mapped_start + 1
required_signal= read_signal[start_idx: end_idx]

For reads mapped to - strand:
start_idx = neg_mapped_end - y- 1
end_idx = neg_mapped_end - x
required_signal= read_signal[start_idx: end_idx]

Note that x and y are inclusive and the length of the required signal = y - x + 1.


text:
Comparing common genes and transcripts between Gencode and Refseq
This is a very naive question - I am trying to compare and get the common lncRNA genes and transcripts between Gencode and Refseq from their gff files.
Since their gene_id, gene_names are different, is there a metric that I can use to compare to see if they share similar genes?
Or a conversion table that would convert one gene/transcript name format into another...
Thank you so much.
gencode.gff3
chr1    HAVANA  gene    29554   31109   .   +   .   ID=ENSG00000243485.5;gene_id=ENSG00000243485.5;gene_type=lncRNA;gene_name=MIR1302-2HG;level=2;hgnc_id=HGNC:52482;tag=ncRNA_host;havana_gene=OTTHUMG00000000959.2
chr1    HAVANA  transcript  29554   31097   .   +   .   ID=ENST00000473358.1;Parent=ENSG00000243485.5;gene_id=ENSG00000243485.5;transcript_id=ENST00000473358.1;gene_type=lncRNA;gene_name=MIR1302-2HG;transcript_type=lncRNA;transcript_name=MIR1302-2HG-202;level=2;transcript_support_level=5;hgnc_id=HGNC:52482;tag=not_best_in_genome_evidence,dotter_confirmed,basic,Ensembl_canonical;havana_gene=OTTHUMG00000000959.2;havana_transcript=OTTHUMT00000002840.1
chr1    HAVANA  exon    29554   30039   .   +   .   ID=exon:ENST00000473358.1:1;Parent=ENST00000473358.1;gene_id=ENSG00000243485.5;transcript_id=ENST00000473358.1;gene_type=lncRNA;gene_name=MIR1302-2HG;transcript_type=lncRNA;transcript_name=MIR1302-2HG-202;exon_number=1;exon_id=ENSE00001947070.1;level=2;transcript_support_level=5;hgnc_id=HGNC:52482;tag=not_best_in_genome_evidence,dotter_confirmed,basic,Ensembl_canonical;havana_gene=OTTHUMG00000000959.2;havana_transcript=OTTHUMT00000002840.1
chr1    HAVANA  exon    30564   30667   .   +   .   ID=exon:ENST00000473358.1:2;Parent=ENST00000473358.1;gene_id=ENSG00000243485.5;transcript_id=ENST00000473358.1;gene_type=lncRNA;gene_name=MIR1302-2HG;transcript_type=lncRNA;transcript_name=MIR1302-2HG-202;exon_number=2;exon_id=ENSE00001922571.1;level=2;transcript_support_level=5;hgnc_id=HGNC:52482;tag=not_best_in_genome_evidence,dotter_confirmed,basic,Ensembl_canonical;havana_gene=OTTHUMG00000000959.2;havana_transcript=OTTHUMT00000002840.1
chr1    HAVANA  exon    30976   31097   .   +   .   ID=exon:ENST00000473358.1:3;Parent=ENST00000473358.1;gene_id=ENSG00000243485.5;transcript_id=ENST00000473358.1;gene_type=lncRNA;gene_name=MIR1302-2HG;transcript_type=lncRNA;transcript_name=MIR1302-2HG-202;exon_number=3;exon_id=ENSE00001827679.1;level=2;transcript_support_level=5;hgnc_id=HGNC:52482;tag=not_best_in_genome_evidence,dotter_confirmed,basic,Ensembl_canonical;havana_gene=OTTHUMG00000000959.2;havana_transcript=OTTHUMT00000002840.1
chr1    HAVANA  transcript  30267   31109   .   +   .   ID=ENST00000469289.1;Parent=ENSG00000243485.5;gene_id=ENSG00000243485.5;transcript_id=ENST00000469289.1;gene_type=lncRNA;gene_name=MIR1302-2HG;transcript_type=lncRNA;transcript_name=MIR1302-2HG-201;level=2;transcript_support_level=5;hgnc_id=HGNC:52482;tag=not_best_in_genome_evidence,basic;havana_gene=OTTHUMG00000000959.2;havana_transcript=OTTHUMT00000002841.2
chr1    HAVANA  exon    30267   30667   .   +   .   ID=exon:ENST00000469289.1:1;Parent=ENST00000469289.1;gene_id=ENSG00000243485.5;transcript_id=ENST00000469289.1;gene_type=lncRNA;gene_name=MIR1302-2HG;transcript_type=lncRNA;transcript_name=MIR1302-2HG-201;exon_number=1;exon_id=ENSE00001841699.1;level=2;transcript_support_level=5;hgnc_id=HGNC:52482;tag=not_best_in_genome_evidence,basic;havana_gene=OTTHUMG00000000959.2;havana_transcript=OTTHUMT00000002841.2
chr1    HAVANA  exon    30976   31109   .   +   .   ID=exon:ENST00000469289.1:2;Parent=ENST00000469289.1;gene_id=ENSG00000243485.5;transcript_id=ENST00000469289.1;gene_type=lncRNA;gene_name=MIR1302-2HG;transcript_type=lncRNA;transcript_name=MIR1302-2HG-201;exon_number=2;exon_id=ENSE00001890064.1;level=2;transcript_support_level=5;hgnc_id=HGNC:52482;tag=not_best_in_genome_evidence,basic;havana_gene=OTTHUMG00000000959.2;havana_transcript=OTTHUMT00000002841.2

refseq.gff
chr1    RefSeq  region  1       248956422       .       +       .       ID=NC_000001.11:1..248956422;Dbxref=taxon:9606;Name=1;chromosome=1;gbkey=Src;genome=chromosome;mol_type=genomic DNA
chr1    Gnomon  lnc_RNA 29774   35418   .   +   .   ID=rna-XR_007065314.1;Parent=gene-MIR1302-2HG;Dbxref=GeneID:107985730,Genbank:XR_007065314.1,HGNC:HGNC:52482;Name=XR_007065314.1;gbkey=ncRNA;gene=MIR1302-2HG;model_evidence=Supporting evidence includes similarity to: 1 mRNA%2C 5 ESTs%2C 10 long SRA reads%2C and 96%25 coverage of the annotated genomic feature by RNAseq alignments%2C including 4 samples with support for all annotated introns;product=MIR1302-2 host gene;transcript_id=XR_007065314.1
chr1    BestRefSeq  lnc_RNA 34611   36081   .   -   .   ID=rna-NR_026818.1;Parent=gene-FAM138A;Dbxref=GeneID:645520,Genbank:NR_026818.1,HGNC:HGNC:32334;Name=NR_026818.1;gbkey=ncRNA;gene=FAM138A;product=family with sequence similarity 138 member A;transcript_id=NR_026818.1
chr1    Gnomon  lnc_RNA 95222   107727  .   +   .   ID=rna-XR_007065337.1;Parent=gene-LOC124903816;Dbxref=GeneID:124903816,Genbank:XR_007065337.1;Name=XR_007065337.1;experiment=COORDINATES: polyA evidence [ECO:0006239];gbkey=ncRNA;gene=LOC124903816;model_evidence=Supporting evidence includes similarity to: 1 EST%2C 8 long SRA reads%2C and 84%25 coverage of the annotated genomic feature by RNAseq alignments%2C including 2 samples with support for all annotated introns;product=uncharacterized LOC124903816%2C transcript variant X1;transcript_id=XR_007065337.1
chr1    Gnomon  lnc_RNA 95222   107727  .   +   .   ID=rna-XR_007065338.1;Parent=gene-LOC124903816;Dbxref=GeneID:124903816,Genbank:XR_007065338.1;Name=XR_007065338.1;experiment=COORDINATES: polyA evidence [ECO:0006239];gbkey=ncRNA;gene=LOC124903816;model_evidence=Supporting evidence includes similarity to: 1 EST%2C 8 long SRA reads%2C and 84%25 coverage of the annotated genomic feature by RNAseq alignments%2C including 2 samples with support for all annotated introns;product=uncharacterized LOC124903816%2C transcript variant X2;transcript_id=XR_007065338.1
chr1    Gnomon  lnc_RNA 120849  175246  .   -   .   ID=rna-XR_001737579.3;Parent=gene-LOC124900384;Dbxref=GeneID:124900384,Genbank:XR_001737579.3;Name=XR_001737579.3;gbkey=ncRNA;gene=LOC124900384;model_evidence=Supporting evidence includes similarity to: 5 long SRA reads%2C and 78%25 coverage of the annotated genomic feature by RNAseq alignments%2C including 34 samples with support for all annotated introns;product=uncharacterized LOC124900384%2C transcript variant X13;transcript_id=XR_001737579.3
chr1    Gnomon  lnc_RNA 146300  175246  .   -   .   ID=rna-XR_007065328.1;Parent=gene-LOC124900384;Dbxref=GeneID:124900384,Genbank:XR_007065328.1;Name=XR_007065328.1;gbkey=ncRNA;gene=LOC124900384;model_evidence=Supporting evidence includes similarity to: 4 ESTs%2C 24 long SRA reads%2C and 76%25 coverage of the annotated genomic feature by RNAseq alignments%2C including 19 samples with support for all annotated introns;product=uncharacterized LOC124900384%2C transcript variant X9;transcript_id=XR_007065328.1
chr1    Gnomon  lnc_RNA 146301  175246  .   -   .   ID=rna-XR_007065325.1;Parent=gene-LOC124900384;Dbxref=GeneID:124900384,Genbank:XR_007065325.1;Name=XR_007065325.1;gbkey=ncRNA;gene=LOC124900384;model_evidence=Supporting evidence includes similarity to: 4 ESTs%2C 25 long SRA reads%2C and 78%25 coverage of the annotated genomic feature by RNAseq alignments%2C including 14 samples with support for all annotated introns;product=uncharacterized LOC124900384%2C transcript variant X6;transcript_id=XR_007065325.1
chr1    Gnomon  lnc_RNA 146435  175246  .   -   .   ID=rna-XR_007065330.1;Parent=gene-LOC124900384;Dbxref=GeneID:124900384,Genbank:XR_007065330.1;Name=XR_007065330.1;gbkey=ncRNA;gene=LOC124900384;model_evidence=Supporting evidence includes similarity to: 1 EST%2C 25 long SRA reads%2C and 75%25 coverage of the annotated genomic feature by RNAseq alignments%2C including 58 samples with support for all annotated introns;product=uncharacterized LOC124900384%2C transcript variant X11;transcript_id=XR_007065330.1
chr1    Gnomon  lnc_RNA 146470  175246  .   -   .   ID=rna-XR_007065322.1;Parent=gene-LOC124900384;Dbxref=GeneID:124900384,Genbank:XR_007065322.1;Name=XR_007065322.1;gbkey=ncRNA;gene=LOC124900384;model_evidence=Supporting evidence includes similarity to: 4 ESTs%2C 24 long SRA reads%2C and 76%25 coverage of the annotated genomic feature by RNAseq alignments%2C including 43 samples with support for all annotated introns;product=uncharacterized LOC124900384%2C transcript variant X3;transcript_id=XR_007065322.1
chr1    Gnomon  lnc_RNA 149034  175246  .   -   .   ID=rna-XR_007065324.1;Parent=gene-LOC124900384;Dbxref=GeneID:124900384,Genbank:XR_007065324.1;Name=XR_007065324.1;gbkey=ncRNA;gene=LOC124900384;model_evidence=Supporting evidence includes similarity to: 4 ESTs%2C 25 long SRA reads%2C and 78%25 coverage of the annotated genomic feature by RNAseq alignments%2C including 17 samples with support for all annotated introns;product=uncharacterized LOC124900384%2C transcript variant X5;transcript_id=XR_007065324.1
chr1    Gnomon  lnc_RNA 149035  175246  .   -   .   ID=rna-XR_007065331.1;Parent=gene-LOC124900384;Dbxref=GeneID:124900384,Genbank:XR_007065331.1;Name=XR_007065331.1;gbkey=ncRNA;gene=LOC124900384;model_evidence=Supporting evidence includes similarity to: 1 EST%2C 24 long SRA reads%2C and 76%25 coverage of the annotated genomic feature by RNAseq alignments%2C including 38 samples with support for all annotated introns;product=uncharacterized LOC124900384%2C transcript variant X12;transcript_id=XR_007065331.1
chr1    Gnomon  lnc_RNA 149035  175246  .   -   .   ID=rna-XR_007065327.1;Parent=gene-LOC124900384;Dbxref=GeneID:124900384,Genbank:XR_007065327.1;Name=XR_007065327.1;gbkey=ncRNA;gene=LOC124900384;model_evidence=Supporting evidence includes similarity to: 4 ESTs%2C 24 long SRA reads%2C and 76%25 coverage of the annotated genomic feature by RNAseq alignments%2C including 33 samples with support for all annotated introns;product=uncharacterized LOC124900384%2C transcript variant X8;transcript_id=XR_007065327.1
chr1    Gnomon  lnc_RNA 151857  175246  .   -   .   ID=rna-XR_007065323.1;Parent=gene-LOC124900384;Dbxref=GeneID:124900384,Genbank:XR_007065323.1;Name=XR_007065323.1;gbkey=ncRNA;gene=LOC124900384;model_evidence=Supporting evidence includes similarity to: 4 ESTs%2C 25 long SRA reads%2C and 81%25 coverage of the annotated genomic feature by RNAseq alignments%2C including 5 samples with support for all annotated introns;product=uncharacterized LOC124900384%2C transcript variant X4;transcript_id=XR_007065323.1


text:
need bam file for pilon
I just ran an assembly on yeast genomes using Flye and I want to polish those assemblies with Pilon but it requires a sorted BAM file.
How do I make a BAM file of the resulting assembled.fasta?

What reads do you have? Only long read or short read or whatever? You need to align the starting reads to your assembly and then supply that bam, as implied by the user172818 answer.

Run a short-read mapper. For example:
bwa index assembled.fasta
bwa mem -pt16 assembled.fasta read1.fq.gz read2.fq.gz \
  | samtools sort -m4G -@4 -o align.bam -

This is good if the user has Illumina, that isn't clear (though it is implied by using pilon).


text:
Adapter trimming
I am trying to do adapter trimming, alignment and sorting for a range of large scale  paired end fastq files. The code I am using is given below:
#!/bin/bash
# Do Trimming, Alignment according to Config file

#SBATCH --job-name=T_Align_single
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=32
#SBATCH --partition=batch
#SBATCH --account=markay
#SBATCH --time=72:00:00
#SBATCH --mem=16G
#SBATCH --mail-type=TIME_LIMIT_80,BEGIN,END
#SBATCH --mail-user=<email-redacted>

module load bowtie2/2.4.1
module load samtools/1.9

mkdir -p Mapping/

config=${INPUT_FILE}

file1=$(cat $config | awk '$1 ~ /fastq/ {print $1}')
file2=$(cat $config | awk '$1 ~ /fastq/ {print $2}')
file3=$(cat $config | awk '$1 ~ /index=/ {print $1}'| sed 's/ref_index=//')

#adapter trimming
trimfile1=Mapping/"${file1%.fastq}".trim.fastq
trimfile2=Mapping/"${file2%.fastq}".trim.fastq

 java -jar /labs/markay/Aranyak/Adriana_intro_CnT/Cnt_Test/Test-analysis/fastq_processing/Trimmomatic-0.39/trimmomatic-0.39.jar  PE -phred33 $file1 $file2 \
   $trimfile1 Mapping/"${file1%.fastq}".unpaired.fastq \
   $trimfile2 Mapping/"${file2%.fastq}".unpaired.fastq \
   LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 \
   MINLEN:20 ILLUMINACLIP: /labs/markay/Aranyak/Adriana_intro_CnT/adapter.fa:2:30:10 -threads ${SLURM_CPUS_PER_TASK}

# mapping using bowtie2
sortbam=Mapping/"${file1%_R1_001.fastq}".sort.bam
sample_ID=$(cat $config | awk '$1 ~ /fastq/ {print $4}')
bowtie2 -p ${SLURM_CPUS_PER_TASK} --local --very-sensitive-local --no-unal \
   --no-mixed --no-discordant --phred33 -I 10 -X 700 \
   -x $file3 \
   -1 $trimfile1 -2 $trimfile2 --rg-id $sample_ID | samtools view -Sb - | samtools sort -@${SLURM_CPUS_PER_TASK} -o $sortbam

#Index the sorting by coordinate bam file
samtools index ${sortbam}

My original fastq files before adapter trimming looked like this
R3-LiD-M-LK_S16_L002_R2_001.fastq.gz
R3-LiD-M-LK_S16_L002_R1_001.fastq.gz

However after adapter trimming I removed the lane names and my file looks like this
R3-LiD-M-LK_S16_R2_001.fastq.gz
R3-LiD-M-LK_S16_R1_001.fastq.gz

The code that worked for the original file is not generating a proper bam file at the end though I am getting a list of trimmed fastq files. The error message I am getting is:
[E::hts_open_format] Failed to open file R1-LiD-H-DJ_S3.sort.rmdup.bam
samtools index: failed to open "R1-LiD-H-DJ_S3.sort.rmdup.bam": No such file or directory
[E::hts_open_format] Failed to open file R1-LiD-H-DJ_S3.sort.bam
The file 'R1-LiD-H-DJ_S3.sort.bam' does not exist
[E::hts_open_format] Failed to open file R1-LiD-H-DJ_S3.sort.rmdup.bam
The file 'R1-LiD-H-DJ_S3.sort.rmdup.bam' does not exist
R1-LiD-H-DJ_S3 DONE.

How should I change the code to generate proper alignment and bam files?

This should be a simple learning exercise for you - all that's changed is that the lane number part is gone, but your error doesn't seem related to that at all, and no variables use the lane number part anywhere. Add a set -x after the #SBATCH lines and look at the job's STDERR file to see the exact commands being executed.

Or just print out your variables so you can see what is happening. You are trying to open a file named R1-LiD-H-DJ_S3.sort.bam and a file named R1-LiD-H-DJ_S3.sort.rmdup.bam and neither of them exist. You also don't mention rmdub anywhere in your code, is this really the code you are running?

Thanks for your help, I have been able to trouble shoot the error with the naming of the files.

Hi @AranyakGoswami, can you please write up your solution as an answer to this question (and make it an accepted answer)? That will help the StackExchange bots do the right thing with regards to not hunting for more answers to an already-solved question.


text:
WCGNA - Relate modules with Y features when the % of variance explained of each eigengen is low
I'm doing a WCGNA analysis (signed network) on microbiome 16S data.
I have transformed counts to centeres log-ratio transformed data (CLR) to address the compositional characteristics of the data and have obtained a pretty decent cluster (around 7 modules, 20 taxa each).
I correlated the different modules eigengenes to my features of interest using both bicor and pearson correlation measures. Overall, the results I'm getting make sense but I'm having some questions regarding the method: The % of variance explained of each eigengen is around 30%. Is this value normal?
Is it ok that I'm using only the PC1 of the module as a summary of the module itself when variances explained of the PC1 are so low?
Thx!!

Which bacterial groups, or perhaps rather what was the original sample for the meta genomics study?

Part 1 Just to clarify on the specific issue of using PC1 alone given the reported result ...
Without knowing the underlying ecology, which is important, 30% is too low to any draw conclusions in PCA-based analysis, particularly because this could be a really complex niche, ecotope, ecoregion, etc ...
Thus this specific result from the unsupervised learning shouldn't be used as the basis for any downstream analysis. This is because the remaining 70% of the variance could overturn PC1, where PC1 could be separating a minor feature of the community ecology.
Thus, its informative because it states the first component is insufficient to proceed - however I stress its difficult to give a clear answer without knowing the expectedly complexity of the meta-genomic sample.
I suspect its a nice result, but in broad general terms its not a good idea.

Part 2

of each eigengen is around 30%.

If PC1-PC3 is 90%, then thats good and maybe its not so complex. Personally, I'd be delighted with that result. However, the inference is that the clustering is not clear. I would simply suggest further data mining is required.

Proportion of variance explained of about 30% is indeed on the low side, but nothing out of the ordinary. I have no experience with 16S data; in gene expression data, between 30 and 50% is actually considered good. PVE byt the eigengene depends strongly on the number of samples; the more samples, the lower the PVE will be.
You can try to calculate additional singular vectors and their PVE. Usually this is much lower than the eigengene (1st singular vector). So I would focus more on the biological significance - do the modules make biological sense?


text:
Recode bgen files against reference genome
I have a set of .bgen files split by chromosomes, including both SNPs and INDELS (named with both rsids and chr_pos_a1_a2). I know a large portion of them are encoded on the wrong strand, and I want to recode them against a reference genome (in my case, 1KGP).
How can I recode the variants in my starting files against the reference genome?

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


text:
Parsing a genbank file and outputting specific feature information to a csv using BioPython
So I am trying to parse through a genbank file, extract particular feature information and output that information to a csv file. The example genbank file looks like this:
SBxxxxxx.LargeContigs.gbk

LOCUS       scaffold_31            38809 bp    DNA              UNK 01-JAN-1980
DEFINITION  scaffold_31.
ACCESSION   scaffold_31
VERSION     scaffold_31
KEYWORDS    .
SOURCE      .
      ORGANISM  .
COMMENT     ##antiSMASH-Data-START##
            Version      :: 6.1.1
            Run date     :: 2022-09-21 11:09:55
            ##antiSMASH-Data-END##
FEATURES             Location/Qualifiers
            protocluster    26198..38809
                            /aStool="rule-based-clusters"
                            /category="terpene"
                            /contig_edge="True"
                            /core_location="[36197:37079](-)"
                            /cutoff="20000"
                            /detection_rule="(Terpene_synth or Terpene_synth_C or
                            phytoene_synt or Lycopene_cycl or terpene_cyclase or NapT7
                            or fung_ggpps or fung_ggpps2 or trichodiene_synth or TRI5)"
                            /neighbourhood="10000"
                            /product="terpene"
                            /protocluster_number="1"
                            /tool="antismash"

Now for the output file, I want to create a csv with 3 columns. One column will have the Scaffold information (ie. scaffold_31), the second column will have the category value in the protocluster feature (ie. /category = "terpene") and the third column will have the product value in the protocluster feature (ie. /product="terpene")
This is what I have so far for code. I know I can sort through the feature.qualifiers in the protocluster feature to get the category and product. I am not sure how to extract the scaffold information. I am completely new to parsing through gene bank files so have little knowledge in this domain. Thanks in advance for any assitance!

import Bio
from Bio import SeqIO
import os

input = "/Path to SBxxxxxx.LargeContigs.gbk"
output = open("output.csv", "w")

if not os.path.exists(output):
     for record in SeqIO.parse(input, "genbank")
          for feature in record.features:
              if feature.type == "protocluster" and "category" and "product" in feature.qualifiers:
                  outfile = feature.qualifiers["category"][0] + "," + feature.qualifiers["product"][0] + "\n"
                  output.write(outfile)

```

Depending on which field you want to pull the "scaffold_31" text from, you have a few options:

GenBank Text
BioPython object

LOCUS
record.name

DEFINITION
record.description

ACCESSION
record.accessions[0]

VERSION
record.version

Python's built in dir() function is handy for figuring out this kind of thing.
If you're working with a draft flat file (like BankIt gives you just before submitting) note that some of those are placeholders that get updated with the actual accession info when it's finalized.  So your "scaffold_31" text will only show up I think in the DEFINITION line in the end if I remember right.

I've used SARS-CoV-2 (Genbank: PA544053), because there was no Genbank entry given in the OPs question. If my example is representative (might not be) I think its about the object attributes. To understand the object I listed its attributes,
print(record[0].__dict__.keys())

Output

dict_keys(['_seq', 'id', 'name', 'description', 'dbxrefs',
'annotations', '_per_letter_annotations', 'features'])

The docs and @jesse's very kind response says there's a 'accession' attribute (Biopython docs below). @Jesse did mention dir() which was cool. I couldn't find record[0].accession or perhaps record[0].accessions and the OP might have had the same problem.
In my example there is an 'annotations' attribute  and beneath that was 'accession' accessed via,
print(record[0].annotations['accessions'])

Output

['PA544053']

There it is ... Here's the full code including the CSV package, I'm using efetch so it'll just copy and paste and run.
from Bio import Entrez
from Bio import SeqIO
import csv
from pathlib import Path

fileout = Path('/mypath/mydir', 'test.csv')
handle = Entrez.efetch(db="nucleotide", id="PA544053", rettype="gb",retmode="text") #,
record = [rec for rec in SeqIO.parse(handle,'genbank')]
fieldnames = ['Accession', 'db_xref', 'organism']

# print(record[0].__dict__.keys())
print(record[0].annotations['accessions'])

with open (fileout, 'w') as f:
    csvwriter = csv.writer(f, delimiter="\t", quotechar='"')
    csvwriter.writerow(fieldnames)
    for feature in record[0].features:
        if feature.type == 'organism' or "mol_type" and "db_xref" in feature.qualifiers:
            csvwriter.writerow([''.join(record[0].annotations['accessions']), ''.join(feature.qualifiers['db_xref']), ''.join(feature.qualifiers['organism'])])

Output:
Accession      db_xref       organism
PA544053       taxon:2697049 Severe acute respiratory syndrome coronavirus 2

Biopython docs
A likely reason for the question is the missing attribute is described in the official docs. Using a GenBank object (not SeqIO) there is certainly an accession attribute, https://biopython.org/docs/1.75/api/Bio.GenBank.html
The official Biopython page ...

from Bio import GenBank
with open("GenBank/NC_000932.gb") as handle:
    record = GenBank.read(handle)
    print(record.accession)

I think the basis of the question is to associate the accession number with the biochemical/genetic info. Note, I don't know the difference between SeqIO and GenBank objects.


text:
convert GT:PL format to GT:GP format
I have a vcf info formatted in GT:PL.
chr22  49994037    .   G   A   345.64  PASS    AC=1;AF=0.500;AN=2;BaseQRankSum=-3.605;DP=35;ExcessHet=3.0103;FS=1.377;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=10.17;ReadPosRankSum=-0.035;SOR=0.392    GT:AD:DP:GQ:PL  0/1:20,14:34:99:353,0,588
I want to convert it into GT:GP. But I have to write a script for that. Is there any way to do it fast and simple?
Note:
GT:Genotype
PL: Phred-scaled Genotype Likelihoods
GP: Genotype Probabilities
check the link below for all the info tags
https://samtools.github.io/hts-specs/VCFv4.2.pdf

Seems like you may be confusing genotype probability and genotype likelihood? That's how you convert GL to PL, not GP to PL.

I think I am right. GP limits to 0 to 1 probability. GL is genotype likelihoods comprised of comma-separated floating point log10-scaled likelihoods and PL is a Phred-scaled likelihood. So, the conversion formulae should be like this. Can u clarify what I have done wrong?

It says PL : the phred-scaled genotype likelihoods, not genotype probabilities. You said you want genotype probabilities, not genotype likelihoods. Why are you trying to convert PL to GP anyway? Genotype probabilities are usually posterior probabilities from imputation and genotype/phred-scaled likelihoods are what you get from variant callers.

to compare my imputation accuracy with true genotype.

So the line you shared is presumably the true genotype, since it doesn't look like it's been imputed? Either way, you want to use dosages for that, not genotype probabilities. Why do you want genotype probabilities?

If u want to use IQS https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0009697#pone-0009697-t002 then u will need all the genotype probabilities. the r squared relation has some disadvantage in lower minor allele frequency

Let us continue this discussion in chat.

Use this tool from GATK to obtain the posterior probabilities and then use bcftools to keep the only the FORMAT fields you want (e.g. GT:GP, ... )
https://gatk.broadinstitute.org/hc/en-us/articles/360036730531-CalculateGenotypePosteriors

Currently you have phred-scaled genotype likelihoods or PL values. Genotype likelihoods give the likelihood of observing the read data (i.e. the number of reads which map to either the ref or alt allele), given a particular underlying genotype (0/0, 0/1, 1/1 at a bi-allelic locus); $GL = p(G|D)$. Phred scaling simply scales it like $PL = -10 * log(GL)$. These are usually given by a genotype caller.
Genotype probabilities are slightly different in that they are posterior probabilities that the genotype is correct, given either other genotypes from a population or trio and a prior distribution for the allele frequencies. Therefore, you can't trivially go from likelihoods to probabilities without accounting for this external information and the prior. As suggested by JRodrigoF, you should use something like CalculateGenotypePosteriors to do this, maybe including 1000 genomes allele counts.


text:
I am trying to use slingshot for pseudotime analysis. I am getting error while saving an object as dataframe
When I run following code, I get error at -
slingshot_df <- data.frame(colData(sce))

the error is :

Error in as.vector(x) : no method for coercing this S4 class to a
vector

library(Polychrome)
library(ggbeeswarm)
library(ggthemes)

# this define the cluster color. You can change it with different color scheme.
my_color <- createPalette(length(levels(sce$ident)), c("#010101", "#ff0000"), M=1000)
names(my_color) <- unique(as.character(sce$ident))

slingshot_df <- data.frame(colData(sce))

# re-order y-axis for better figure: This should be tailored with your own cluster names
# slingshot_df$ident = factor(slingshot_df$ident, levels=c(4,2,1,0,3,5,6))

ggplot(slingshot_df, aes(x = slingPseudotime_1, y = ident,
                              colour = ident)) +
    geom_quasirandom(groupOnX = FALSE) + theme_classic() +
    xlab("First Slingshot pseudotime") + ylab("cell type") +
    ggtitle("Cells ordered by Slingshot pseudotime")+scale_colour_manual(values = my_color)


text:
Visualize RNA co-folded secondary structure with Python
I would like to visualize the secondary structure of co-folded RNA strands in Python, preferably in a Jupyter notebook.
What are the recommended tools? Can anyone provide an example?

... . I know both aspects of this question and I'm not aware of an interactive Python package using e.g. iPython magic commands that will do this. The thing you could try is  draw_RNA https://github.com/DasLab/draw_rna
draw_RNA will operate on Jupyter, but I've not tried it.
A very, very good package is VARNA (its Java BTW). The VARNA group make further recommendations from there for complex structures. I would vouch for this package for complex structures.
A good approach is to output in Stockholm 1.0 format. It then just drops and drags into VARNA. Dot-bracket format also works and I would assume ct format works too.
A good reason why in-line graph imaging isn't the best solution for complex structures, i.e. the in-line graphs in a Python debugger or Jupiter is complex RNA structures need to be interactive because you need to move them: they "clump" together and have to manually separate them.
As a final note Forna may be of benefit,


text:
Cannot run `loadData()` after saving RDS file and reading RDS file using `saveRDS` & `readRDS` in Seurat
I want to load multiple sample data(pbmc16 and pbmc45), and integrate by following the tutorial below

Performing integration on datasets normalized with SCTransform
https://satijalab.org/seurat/articles/integration_introduction.html
performing integration on datasets normalized with sctransform-1

What I did:
1. I performed the pre-processing step (filtering) and SCtransform in each sample, and then save each data with the code
saveRDS(pbmc16, file = "../output/pbmc16_Seurat_Guided_Clustering_Tutorial.rds")
saveRDS(pbmc45, file = "../output/pbmc45_Seurat_Guided_Clustering_Tutorial.rds")
> pbmc16; type S4 [20542x20700](SeuratObject:: Seurat) assays; list [2]
> meta.data; list [20700x8] (S3:data.frame)

2. read each sample data by using readRDS function
pbmc16 <- readRDS("../output/pbmc16_Seurat_Guided_Clustering_Tutorial.rds")
pbmc45 <- readRDS("../output/pbmc45_Seurat_Guided_Clustering_Tutorial.rds")

3. After that, I want to load the data as shown in the tutorial
(LoadData("ifnb"))
LoadData(pbmc16)
LoadData(pbmc45)

Then I got this error
> LoadData(pbmc16) Error in match(x, table, nomatch = 0L) : 'match'
> requires vector arguments

Problem: I don't know how to load each data and integrate the multiple data.
Could you please suggest me the recommended code or advice?

Hi, you seem to be following the tutorial on the web which unfortunately is not up to date with the package code (in current manual in the web doesn't seem to exists a LoadData function).


text:
KEGG retrieve hierarchy for organism in biopython
I want to obtain a direct graph from KEGG as well as the genes present in a single term. While it is pretty easy to get list of genes/compounds that are part of a specific category, I am kinda stucked to obtain the hierarchical relationships between terms. The docs on biopython.KGML are not really helpful. This is what I have so far
from Bio import SeqIO
import Bio.KEGG.REST as kegg
from Bio.KEGG.KGML import KGML_parser
from Bio.Graphics.KGML_vis import KGMLCanvas
from Bio.Graphics.ColorSpiral import ColorSpiral
from itertools import chain

k =  kegg.kegg_list('pathway', 'pae').read()
todf = []
for x in k.split('\n'):
    if x=='':
        next
    dd = x.split(' - ')
    dd = [v.split('\t') if '\t' in v else v for v in dd]
    dd = list(chain.from_iterable(dd))[:2]
    try:
        ll= kegg.kegg_get(dd[0]).read().split('GENE')[-1]
        ll = [x.lstrip() for x in ll.split('COMPOUND')[0].split('\n')]
        ll= [x.rstrip().split('\t') for x in ll]
        for gn in ll:
            gn.extend(dd)
        todf.append(pd.DataFrame(ll))
    except IndexError:
        pass
todf = pd.concat(todf)
# remove compounds and random strings
todf= todf[todf[0].str.startswith('PA')]
todf[0] = [x.split(';')[0] if ';' in x else x for x in todf[0]]
todf[0] = todf[0].str.split('  ')
todf=todf.explode(0)

and now I have the protein to term relationship but not sure how to get the hierarchy between terms.
Any idea?

Ok realized it was much easier than expected.
On KEGG it is possible to DL a json file, which is intrinsically ordered. So   this parses it into a df
import pandas as pd
import json

todf = []

kegg = json.load(open('meta/KEGG_pae00001.json'))
for main in kegg['children']:
    for broad in main['children']:
        for sub in broad['children']:
            for gn in sub.get('children', [None]):
                if gn:
                     nm = gn['name'].split(';')[0].split(' ')
                     todf.append([nm[0], nm[1], sub['name'], broad['name'], main['name']])
                else:
                     # some terms have no children
                     pass
                     #print(sub)

todf = pd.DataFrame(todf)
todf.columns=['Locus', 'GeneName', 'Cat1', 'Cat2', 'Cat3']
todf.to_csv('meta/kegg_parsed.csv')

Hi, assuming this answer works for you, could you please mark this as the accepted answer? That will help the automatic answer-finding algorithms do the right thing.


text:
How can I assemble my genome from raw files?
I've had my whole genome sequenced (at 30x average coverage) by a lab, and they have provided the raw files to me (BAM, FASTQ, and VCF).
How can I assemble it?
And does assembly provide any further insights that I can't get from the separated fragments?
I'd like to assemble and annotate it, and find insights, but I don't know where to start.

What do you want to achieve by 'assembling' it?  That's usually restricted to random organisms where we don't have a reference sequence. Your VCF/BAM/FASTQ data has already been mapped to an assembled genome so it doesn't really make sense to do an assembly of it. You could certainly annotate it thought, but thats a different question altogether.

It would be useful to know the reference genome.

Given you've tagged this as 'human-genome' there is likely little value in performing genome assembly. De-novo assembly of WGS is generally for non-model organisms without a reference genome available, and would require much higher than 30X coverage (and ideally long reads). By providing the BAM and VCF, they'e already saved you some upstream work unless there are needs to re-process the data from FASTQ

No, assembling won't tell you anything new. When working with organisms whose reference genome is known, all you care about is where your sequences differ from the reference. If you know where your sequences differ, then you can infer the rest, they will be whatever is in the reference.
All you want to do is identify variants, those sites where your genome differs from the reference. That's what the VCF file contains. Bear in mind that the reference isn't some sort of platonic ideal of the perfect human, it's just the collected, pooled DNA of 13 anonymous, apparently healthy individuals. You can expect anywhere from about 3 to 5 million differences in your own genome, the overwhelming majority of which will be completely innocuous.
Since you already have your VCF file, your list of variants, that's all you need to investigate. If you really want to, you can make your assembly by taking the reference genome and then adding the variants in your VCF file, but it won't be of much practical use.

To add the "do not assemble your genome" answers. There might be specific circumstances when you would actually like to do that. For example, if the 30x you got are HiFi long reads and have a reason to suspect there are some interesting large rearrangements in your genome (rare disease?), it might be informative to actually assemble the genome and compare it to the reference to find these rearrangements.
Of course there is also plenty of possibilities in between (still assuming long reads), for example to map your sequencing data to a pan-genome graph or simply use a mapper that is aware of genomic rearrangements (like ngmlr) and call the structural variants.
But in the end it boils down to two aspects you did not mention in your question: What would you like to learn about your genome?


text:
Searching for HLA-B in DNA results
I'm trying to find the HLA-B*15:01 variant in my DNA results, prompted from this research paper: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8142661/#:~:text=HLA%2DB*15%3A01%20is%20strongly%20associated%20with%20asymptomatic,mechanism%20underlying%20early%20viral%20clearance.
I have DNA test results from Dante Labs that contain the following files
"*.bam"
"*.bam.bai"
"*.cnv.vcf.gz"
"*.cnv.vcf.gz.tbi"
"*.filtered.indel.vcf.gz"
"*.filtered.indel.vcf.gz.tbi"
"*.filtered.snp.vcf.gz"
"*.filtered.snp.vcf.gz.tbi"
"*1_SA_L001_R1_001.fastq.gz"
"*SA_L001_R2_001.fastq.gz"

I tried loading this on the Genome Explorer on Sequencing.com, but it couldn't find HLA-B at all, only HLA-A and HLA-C.
Then I tried loading the VCF files in IGV. I was able to filter to Chromosome 6 and HLA-B, which resolved to position
chr6:31,319,652-31,326,956

The CNV VCF says no variables found for chr6:31,319,652-31,326,956
The filtered Indel VCF did produce some variants for this position, but I have no clue how to search deeper specifically for HLA-B*15:01. Here is a screenshot:
The filtered SNP VCF produced alot more variants for this position, but again, I don't know how to go deeper than this.
I'm not really sure what to do from here?
I don't have data on what the reference was, so that is unknown to me. I will reach out to Dante Labs about that, is there some sort of standard that is supposed to be used, or is it different for every genome provider?
How do I actually perform the check in IGV assuming the reference did not have it and the VCF did. I also assume this VCF does not have this allele, but I frankly am lost into how to confirm it to the "last mile" through the VCF so to speak, or even just go straight to BAM

Typical SNP calling wouldn't give you good results on classical HLA genes. Ignore the VCF. Try a dedicated HLA genotyper on raw reads instead. I heard optitype is a good one for exome or whole-genome data. I haven't used it myself, though.

The specific question and the easiest question is the absence of HLA-B*15:01 from your sample and that is not unexpected at all.
The maximum frequency of this allele is 40% found in "China Beijing pop 2" , this however is really unusually high and was published by Guang Yang et al (2009, IJI 37, 205-212).
Mostly this allele is present at 5% for moderate frequency of occurrence for example Europe, such as in Germany here. Moreover, in some populations the frequency is very low <<1%, e.g. India or Israel.
So a sample genome at random the will not have the B*15:01.

The more complicated issue is the VCF and the reference genome. This is because the basis on whether an allele is present or absent is a comparison against the reference genome which will generate the VCF. Thus understanding the HLA-B alleles of the reference genome is really important. You need a reference genome which lacks HLA-B*15:01, therefore if your VCF also lacks the allele, its missing. However, if HLA-B*15:01 is present in the reference genome your VCF will not present the allele if its present.
The probability is that the allele is missing from both the reference and sample genome.

Comments
Easy way The easiest way is search the IGV for LT618832.1 (HLA-B*15:01 reference), IGV has a motif finder, here https://software.broadinstitute.org/software/igv/motif_finder
Hard way
This is how I'd do it ..

convert the bam to fasta then
run a blast against Genbank reference LT618832.1

I'd section the bam file into the locus of interest (e.g. chromosome 6) e.g. ... samtools view ... ChrName: ... > chromo6.bam: you'd need to cross check how to do this last bit ... then via samtools >1.3,
samtools fasta chromosome6.bam > chromosome6.fasta

Finally you download a standalone version of blast, I assume its on -c bioconda ...
blastn -query HLAB1501.faa -db chromosome6.fa -out myhits.txt

From memory I don't think you need to makeblastdb for a fasta onto a fasta db.
You'd need to map the SNPs in introns parsing (the coordinates are in the Genbank reference), or optically against the Genbank reference onto the blast output.
You may want to cross check the HLA-B*15:01 against https://www.ncbi.nlm.nih.gov/nuccore/?term=HLA-B%2A15%3A01
[but be careful about the introns]
motif_finder seems easier and accommodates potential heterozygosity.


text:
Writing to .cif from BioPandas
My question is twofold. First, I've run into an issue working with .cif files and pandas dataframes. With BioPandas, I can read a .pdb file into a dataframe and then go the opposite direction, writing a .pdb file from a dataframe. With .cif files, I can only write into a dataframe - PandasMmcif does not contain a to_cif function like PandasPdb has a to_pdb function. Are there any alternative methods I could use to get from .cif to pandas dataframe back to .cif?
Second, the reason I am not simply converting all of my .cif files to .pdb is that I am uncertain whether I will lose data in the conversion. Is this fear unfounded? I don't know the specific differences between .pdb and .cif.

You may lose some metadata during mmCIF <-> PDB conversion (what is converted and what is ignored depends on the converter that you use), but it's the same when you load a PDB or mmCIF file into a program, modify it and write it back – some metadata may get lost. The most essential things should be preserved. Unless you have a big mmcif file that can't be converted to a standard PDB due to the limits on the number of chains and atoms in the PDB format (there are workarounds for this used by some converters).

Perfect, thank you!


text:
Nested variable expansion in WDL code block
I am developing a gatk pipeline for variant calling and want to include trimmomatic as a preprocessing step. I would like the user to specify a program file of parameters to feed into trimmomatic, and a file of adapter sequences to provide to the ILLUMINACLIP parameter. My goal is to include a variable (“~{adapter_file}” below) in the program file that will expand to include the adapter filename, like so:
Adapter file:
>IL_Multiplex_five_prime_fwd
AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT
>IL_Multiplex_three_prime_rev
AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC

Program file:
ILLUMINACLIP:~{adapter_file}:2:30:10:7:keepBothReads
MINLEN:36
SLIDINGWINDOW:4:15
LEADING:3
TRAILING:3

In this way, the user can change the adaptor file without having to alter the program file each time. However, I can’t seem to get the ~{adapter_file} variable to expand inside the ~{trim_program_array} variable. Here’s a simplified version of my trimmomatic WDL file:
Sample workflow:
version 1.0

Workflow trim {
    input {
        File fastq1 = “R1.fastq.gz”
        File fastq2 = “R2.fastq.gz”
        File program_file = “trim_program.txt”
        File adapter_file = “adapters.fa”
    }

    Array[String] trim_program_array = read_lines(program_file)

    call trimmomatic {
        input:
            fastq1 = fastq1,
            fastq2 = fastq2,
            trim_program_array = trim_program_array,
            adapter_file = adaptor_file
    }
}

Task trimmomatic {
    input {
        File fastq1
        File fastq2
        File outfile1 = “trimmed_R1.fastq.gz”
        File outfile2 = “trimmed_R2.fastq.gz”
        Array[String] trim_program_array
        File adapter_file
     }

     command {
        TrimmomaticPE \
        ~{fastq1} ~{fastq2} \
        ~{outfile1} ~{outfile2} \
        ~{sep = “ “ trim_program_array}
     }

     output {
        File out1 = outfile1
        File out2 = outfile2
     }
}

My goal is for the first element of ~{sep = " " trim_program_array} (aka the first line of trim_program.txt) to evaluate to “ILLUMINACLIP:adapters.fa:2:30:10:7:keepBothReads” rather than “ILLUMINACLIP:~{adapter_file}:2:30:10:7:keepBothReads”. I know this is because ~{adapter_file} is nested inside trim_program_array[0], so it’s only expanding the “outer” variable. How can I convince WDL to expand the variable nested inside the first element of that array before executing the bash command?

Hi, I've corrected your question to make the spelling of the word 'adapter' consistent, but noticed that some of those changes were in the configuration file. Can you please re-check to make sure this spelling inconsistency is not causing your problems?

The ~{} syntax only works within the wdl command section or wdl string literals. Since the string is getting read from a file in this case, it is probably getting escaped. It looks like this file is just getting used to provide arguments in the command. Can't you just include the contents of the program file in the command section directly?

@DavyCats Yes, that's what I eventually figured out, but it would be nice to have it in the program file so that users of the pipeline aren't left wondering what ":2:30:10:7:keepBothReads" in the first line of the program file refers to, since 'ILLUMINACLIP:~{adaptor_file}' is missing. There are worse problems to have though, so I'm fine with that work-around. Thanks for the suggestion.

@DavyCats Upon thinking about this further, I realized that one advantage to supplying a program file to Trimmomatic is that the user gets to determine which parameters to trim by, rather than just choosing the value for each predetermined parameter. By encoding "ILLUMINACLIP:~{adaptor_file}" into the code block, I'm in effect making that a required parameter. I've seen many cases where a user wouldn't want to supply that parameter to Trimmomatic at all. I'm hoping that someone will come along and help me with a better solution to this problem so that I don't have to hard code any parameters.

Putting the placeholder inside of an input file won't work, since the contents will just be turned into a string when read with read_lines and won't get treated like WDL code. Since the lines from the program_file are simply getting pasted into the command, you can move the contents into the command section. In a placeholder (~{}), adding a string to an optional value will result in an empty string if that optional value is not supplied (see the spec here). So you can use optional inputs to make the different lines non-mandatory:
task trimmomatic {
    input {
        File fastq1
        File fastq2
        String outfile1 = “trimmed_R1.fastq.gz”
        String outfile2 = “trimmed_R2.fastq.gz”
        File? adapter_file
        Int? minlen
        String? slidingwindow
        Int? leading
        Int? trailing
     }

     command {
        TrimmomaticPE \
        ~{fastq1} ~{fastq2} \
        ~{outfile1} ~{outfile2} \
        ~{"ILLUMINACLIP:" + adapter_file + ":2:30:10:7:keepBothReads"} \
        ~{"MINLEN:" + minlen} \
        ~{"SLIDINGWINDOW:" + slidingwindow} \
        ~{"LEADING:" + leading} \
        ~{"TRAILING:" + trailing}
     }

     output {
        File out1 = outfile1
        File out2 = outfile2
     }
}

If, say, adapter_file were to be omitted in the inputs, then the produced command should look something like this:
TrimmomaticPE \
R1.fq R2.fq \
trimmed_R1.fastq.gz trimmed_R2.fastq.gz \
 \
MINLEN:36 \
SLIDINGWINDOW:4:15 \
LEADING:3 \
TRAILING:3

Also note that the output file paths in the input section should be Strings. Since they don't exist yet, the execution engine will complain that it can't find them if they are Files
Ps. I'm not familiar with the actual trimmomatic command options, so please do check that this produces the desired command. I suspect that the :2:30:10:7:keepBothReads bit could be split up into additional inputs as well.


text:
select highly variable genes out of dataframe
lets say I have the following dataframe:
Ensembl_ID          A   B   C   D   E   F   G
ENSG00000236601.2   0   0   0   1   0   0   1
ENSG00000237094.12  2   187 442 109 88  144 486
ENSG00000269732.1   0   0   0   0   0   0   0
ENSG00000284733.1   0   0   0   0   0   0   0
ENSG00000233653.3   0   0   0   0   0   0   0
ENSG00000250575.1   1   11  17  5   5   8   17
ENSG00000278757.1   0   0   2   0   1   0   7
ENSG00000230021.9   0   25  30  27  9   8   11
ENSG00000235146.2   0   0   0   0   0   0   0
ENSG00000225972.1   0   7   3016    11  5   14  5
ENSG00000225630.1   1   113 194 76  91  47  94
ENSG00000237973.1   25  1037    9767    1160    886 1321    5220

What I would like to do is only keep the ~top 50 highly variable genes:
I thought it would be like this:  dif_genes <- order(apply(df, 1, var), decreasing=TRUE)[1:50]
but that leaves me with a vector containing variances without genes.

One way is adding a new var column to your original data frame and sorting the updated data frame by the newly introduced column:
df$var <- apply(df, 1, var)
df[order(df$var, decreasing=TRUE),]

the OP wants top 50. I would suggest df[order(apply(df[, -1], 1, var), decreasing=TRUE),][1:50, ] to avoid adding the variance as a sample (in case OP wants to do a heatmap) and removing the names that seem to be in the first column. The key @user16197 to keep the top 50 genes is the comma after the 50; it tells R that this is a 2D data.


text:
How to improve this WGCNA analysis
This is my first time working on using WGCNA on a microarray dataset. One of the major problems I am facing is merging close modules which is not really working well. I have quantile normalised the data before working on it. I have put my code below and link to some of the plots generated during this process. I was wondering can somebody advice on what could be improved upon in this analysis?
Link to plots: ![plots] https://i.stack.imgur.com/dF6NA.jpg
I have a total of 767 samples and 13466 probes.
powers = c(c(1:10))
## Scale free topology since it can find hubs
sft = pickSoftThreshold(bryois_norm_keep_use, powerVector = powers, verbose = 5)
par(mfrow = c(1,2));
cex1 = 0.9;
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
labels=powers,cex=cex1,col="red");
abline(h=0.90,col="red")
# Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],
xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
# Turn data expression into topological overlap matrix
power=sft$powerEstimate
# Even though estimate suggests 3, trying ahigh number
# Turn adjacency into topological overlap matrix (TOM)
adjacency <- adjacency(bryois_norm_keep_use, power = 6)
TOMadj <- TOMsimilarity(adjacency,verbose = 5)
dissTOMadj <- 1- TOMadj
# Clustering using TOM
# Call the hierarchical clustering function
hclustGeneTree <- hclust(as.dist(dissTOMadj), method = "average")
# Plot the resulting clustering tree (dendogram)
sizeGrWindow(12, 9)
plot(hclustGeneTree, xlab = "", sub = "",
main = "Gene Clustering on TOM-based disssimilarity",
labels = FALSE, hang = 0.04)
# Make the modules larger, so set the minimum higher
minModuleSize <- 30
# Module ID using dynamic tree cut
dynamicMods <- cutreeDynamic(dendro = hclustGeneTree,
distM = dissTOMadj,
deepSplit = 2, pamRespectsDendro = FALSE,
minClusterSize = minModuleSize,verbose = 5)
table(dynamicMods)
# Convert numeric lables into colors
dynamicColors <- labels2colors(dynamicMods)
table(dynamicColors)
# Plot the dendrogram and colors underneath
sizeGrWindow(8,6)
plotDendroAndColors(hclustGeneTree, dynamicColors, "Dynamic Tree Cut",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05,
main = "Gene dendrogram and module colors")
# Calculate eigengenes
dynamic_MEList <- moduleEigengenes(bryois_norm_keep_use, colors = dynamicColors)
dynamic_MEs <- dynamic_MEList$eigengenes
library(dynamicTreeCut)
# Calculate dissimilarity of module eigengenes
dynamic_MEDiss <- 1-cor(dynamic_MEs)
dynamic_METree <- hclust(as.dist(dynamic_MEDiss))
# Plot the hclust
sizeGrWindow(7,6)
plot(dynamic_METree, main = "Dynamic Clustering of module eigengenes",
xlab = "", sub = "")
######################## MERGE SIMILAR MODULES
dynamic_MEDissThres <- 0.95
# Plot the cut line
#abline(h = dynamic_MEDissThres, col = "red")
# Call an automatic merging function
merge_dynamic_MEDs <- mergeCloseModules(bryois_norm_keep_use, dynamicColors, cutHeight = dynamic_MEDissThres, verbose = 5)
# The Merged Colors
dynamic_mergedColors <- merge_dynamic_MEDs$colors
# Eigen genes of the new merged modules
mergedMEs <- merge_dynamic_MEDs$newMEs
mergedMEs
table(dynamic_mergedColors)
sizeGrWindow(12,9)
plotDendroAndColors(hclustGeneTree, cbind(dynamicColors, dynamic_mergedColors),
c("Dynamic Tree Cut", "Merged dynamic"),
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)

I took a quick look and don't see anything wrong with the module merging. You could lower the merging threshold a bit - 0.95 is very high, I would use 0.9 or 0.85. Other than that though, really nothing seems out of the ordinary.

Hi Peter, Sorry for the late reply but thank you again for your feedback. :)


text:
Make mpileup file with several BAM files
I'm new to this business, so my question may seem as dumb. I want to make an mpileup file using a reference sequence and several (two to begin with) aligned sorted BAM files. And I want the output file to have certain easy-to-read tab-delimited formatting (which is helpful for subsequent analyses of the output file):
field 1 - chromosome,
field 2 - position,
field 3 - reference base,
field 4 - base(s) in SAMPLE1.BAM for position in field 2,
field 5 - number of reads supporting bases in field 4,
field 6 - base(s) in SAMPLE2.BAM for position in field 2,
field 7 - number of reads supporting bases in field 6.
Something like shown below:
field1   field2   field3   field4   field5   field6   field7
Chr1     15967    T        T,A      15,17    T,_      25,_

where _ means no alternative (or second) base (or zero reads) for this position in the BAM file. '0' or empty space (or nothing) instead of _ is ok. It's ok if information in fields 4-7 is shown not in a single field but in several fields.
Which command is suitable for such a purpose?
Thank you.

Have you tried using samtools mpileup in.bam [in2.bam [...]]? Here's information about the output format from the manual page:

The first three columns give the position and reference:

Chromosome name.
1-based position on the chromosome.
Reference base at this position (this will be “N” on all lines if -f/--fasta-ref has not been used).

The remaining columns show the pileup data, and are repeated for each input BAM file specified:

Number of reads covering this position.
Read bases.  This encodes information on matches, mismatches, indels, strand, mapping quality, and starts and ends of reads.

This format seems nearly identical to what you've asked for, so hopefully is close enough for you to tweak your downstream analyses.


text:
Filtering reads greater than 5 from HT-seq count files
I have some raw counts from HTSeq after aligning with hg38 human reference genome. I want to do filtering in a way that the filtered count files should have the same number of lines. The reason behind this is that I want to use DESeq2 for differential gene expression which requires all the count files to be of the same length to carry out the analysis.
I have written a python program to select reads greater than 5 from the second column.
My count files look like this
ENSG00000000003 20
ENSG00000000005 0
ENSG00000000419 123
ENSG00000000457 35
ENSG00000000460 56
ENSG00000000938 0
ENSG00000000971 0
ENSG00000001036 164
ENSG00000001084 83
ENSG00000001167 68
ENSG00000001460 8
ENSG00000001461 13
ENSG00000001497 80
ENSG00000001561 0
ENSG00000001617 39
ENSG00000001626 0
ENSG00000001629 146
ENSG00000001630 4

Python program I have written is as follows
import sys

inputlist = [sys.argv[1], sys.argv[2], sys.argv[3], sys.argv[4], sys.argv[5], sys.argv[6], sys.argv[7], sys.argv[8]]

for i in range(len(inputlist)):
    inputfile = open(inputlist[i], 'r')
    for line in inputfile:
        splitline = line.strip().split('\t')
        if int(splitline[1]) >= 5:
            if not splitline[0] in genelist:
                genelist.append(splitline[0])
    inputfile.close()

outputlist = open(sys.argv[9], 'w')
for gene in genelist:
    outputlist.write(gene + '\n')
outputlist.close()

for i in range(len(inputlist)):
    inputfile = open(inputlist[i], 'r')
    outputfile = open(sys.argv[i+10], 'w')
        for line in inputfile:
        splitline = line.strip().split('\t')
        if splitline[0] in genelist:
            outputfile.write(line)
    inputfile.close()
    outputfile.close()

python filtercounts.py 'mir99-1-L1count.txt' 'mir99-1-L2count.txt' 'mir99-2-L1count.txt' 'mir99-2-L2count.txt' 'mir99-3-L1count.txt' 'mir99-3-L2count.txt' 'genelist.txt' 'mir99-1-L1filtered.txt' 'mir99-1-L2filtered.txt' 'mir99-2-L1filtered.txt'  'mir99-2-L2filtered.txt'  'mir99-3-L1filtered.txt' 'mir99-3-L2filtered.txt'
This program has worked well for me in past , but now I see that although we see many such reads removed but there are many zero reads or reads less than 5 in the dataset.
Some of my output looks like this
ENSG00000000003 0
ENSG00000000419 0
ENSG00000000457 15
ENSG00000000460 12
ENSG00000001036 0
ENSG00000001084 0
ENSG00000001167 0
ENSG00000001460 1
ENSG00000001461 0
ENSG00000001497 1
ENSG00000001617 0
ENSG00000001629 0
ENSG00000001631 0
ENSG00000002016 2
ENSG00000002330 2
ENSG00000002549 0
ENSG00000002822 0
ENSG00000002834 0
ENSG00000002919 0
ENSG00000003056 4
ENSG00000003147 0
ENSG00000003249 1

Any help regarding the problematic behavior of the problem will be appreciated.

After seeing that you want to use DESeq2, I'd recommend you do the count filtering within R, after generating the count matrix. Note that DESeq2 includes functions for loading from HTSeq-count output files:
https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#htseq-count-input
It is not necessary to filter out low count data when processing with DESeq2, but if desired this can be done after converting to a DESeq2 object:

While it is not necessary to pre-filter low count genes before running the DESeq2 functions, there are two reasons which make pre-filtering useful: by removing rows in which there are very few reads, we reduce the memory size of the dds data object, and we increase the speed of the transformation and testing functions within DESeq2. It can also improve visualizations, as features with no information for differential expression are not plotted.

keep <- apply(counts(dds), 1, function(x){all(x >= 5)} # for all gene counts >= 5
keep <- apply(counts(dds), 1, function(x){any(x >= 5)} # for any gene counts >= 5
dds <- dds[keep,]

https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#pre-filtering
Note that DESeq2 cares about replicates, and will almost always pin an adjusted p-value near 1 if any conditions have 1 or 2 replicates (three is the minimum number required to identify outlier values). There's not really any way to work around that other than including more replicates for the low-replicate conditions.
Another program might give different results, but that won't get around the fundamental issue: two replicates is not sufficient to determine sample variation. Either accept that p-values cannot be used as a reliable estimate of model fit this case, or add more replicates.

Should I use any other tool like edgeR which I have not used previously?

@AranyakGoswami, I've added additional clarification about p-values. The replicate numbers simply aren't sufficient to get reliable distribution estimates.


text:
Does cutadapt trim trailing N's first and then use max_n to filter reads?
Background I want to trim leading and (likely just) trailing N's from my WES (Illumina NextSeq500) reads with cutadapt (--trim-n).
I also want to filter for reads with > 40%? N's (--max-n).
Question Does cutadapt first remove the trailing/leading N's before removing reads with >X% N's? This is because if, hypothetically, 40% of the N's are on the end of the reads, I want to keep them.
Cheers and thanks in advance.

Yes, it behaves just like you expect, where a read with >40% Ns but only at the ends will still be kept (after trimming).

Cutadapt runs in two stages, with read modification first and read filtering after.  So first, --trim-n removes Ns at the start and end of reads are removed, and only then is a rule like --max-n applied.
Here's an example where the first 5 sequences have Ns at the ends and the next 5 have Ns interspersed:
>seq01  0%
CACATCAAGAGGTCCTTGGCGACTATTAGCGCCCGGCGGACCCAATCGTA
>seq02 20%
NNNNNCAAGAGGTCCTTGGCGACTATTAGCGCCCGGCGGACCCAANNNNN
>seq03 40%
NNNNNNNNNNGGTCCTTGGCGACTATTAGCGCCCGGCGGANNNNNNNNNN
>seq04 60%
NNNNNNNNNNNNNNNTTGGCGACTATTAGCGCCCGNNNNNNNNNNNNNNN
>seq05 80%
NNNNNNNNNNNNNNNNNNNNGACTATTAGCNNNNNNNNNNNNNNNNNNNN
>seq06 10%
CACATCANGAGGTCCTTNGCGACTATTNGCGCCCGGCNGACCCAATCNTA
>seq07 20%
CANATCANGAGGNCCTTNGCGANTATTNGCGCNCGGCNGACCNAATCNTA
>seq08 40%
CNNATCNNGAGNNCCTNNGCGNNTATNNGCGNNCGGNNGACNNAATNNTA
>seq09 60%
CNNANCNNGNGNNCNTNNGNGNNTNTNNGNGNNCNGNNGNCNNANTNNNA
>seq10 80%
CNNNNNNNGNNNNCNNNNGNNNNTNNNNGNNNNCNNNNGNNNNANNNNNA

With just read filtering:
$ cutadapt --quiet --max-n 0.4 seqs.fa -o -
>seq01  0%
CACATCAAGAGGTCCTTGGCGACTATTAGCGCCCGGCGGACCCAATCGTA
>seq02 20%
NNNNNCAAGAGGTCCTTGGCGACTATTAGCGCCCGGCGGACCCAANNNNN
>seq03 40%
NNNNNNNNNNGGTCCTTGGCGACTATTAGCGCCCGGCGGANNNNNNNNNN
>seq06 10%
CACATCANGAGGTCCTTNGCGACTATTNGCGCCCGGCNGACCCAATCNTA
>seq07 20%
CANATCANGAGGNCCTTNGCGANTATTNGCGCNCGGCNGACCNAATCNTA
>seq08 40%
CNNATCNNGAGNNCCTNNGCGNNTATNNGCGNNCGGNNGACNNAATNNTA

With trimming also:
$ cutadapt --quiet --max-n 0.4 --trim-n seqs.fa -o -
>seq01  0%
CACATCAAGAGGTCCTTGGCGACTATTAGCGCCCGGCGGACCCAATCGTA
>seq02 20%
CAAGAGGTCCTTGGCGACTATTAGCGCCCGGCGGACCCAA
>seq03 40%
GGTCCTTGGCGACTATTAGCGCCCGGCGGA
>seq04 60%
TTGGCGACTATTAGCGCCCG
>seq05 80%
GACTATTAGC
>seq06 10%
CACATCANGAGGTCCTTNGCGACTATTNGCGCCCGGCNGACCCAATCNTA
>seq07 20%
CANATCANGAGGNCCTTNGCGANTATTNGCGCNCGGCNGACCNAATCNTA
>seq08 40%
CNNATCNNGAGNNCCTNNGCGNNTATNNGCGNNCGGNNGACNNAATNNTA


text:
How to use a .out Kraken output format
I was given some Kraken output files but the file extension is .out. I know how to use .KREPORT files but not the .out file. Is there any way convert it to a KREPORT file or what do you recommend?
Thanks in advance.

The extension of a file does not determine its content. What is the format of your .out files?

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.

If I read the file using Notepad++ I can see something like:
C   A00700:50 342 120|120 0:79 821:1 0:36 |:| A:1 0:58 821:4 0:53
U   A00700:50 342 139|106 0:69 226186:4 0:7 226186:5 0:31 |:| A:1 0:25 818:3 0:45
Is it useful?
Thanks.


text:
Identifying highly variable genes in scRNA-seq: Seurat vs M3Drop
Following a single-cell RNA-seq workshop, I created a Seurat object (my_data), normalized the data, and then tried to identify highly variable genes using two different R packages: Seurat and M3Drop.
variable_genes_Seurat <- my_data %>%
  FindVariableFeatures(selection.method = 'vst', nfeatures = 2000) %>%
  VariableFeatures()

variable_genes_M3Drop <- my_data %>%
  GetAssayData('counts') %>% # unnormalized
  NBumiConvertData() %>%
  NBumiFitModel() %>%
  NBumiFeatureSelectionCombinedDrop(ntop = 2000) %>%
  rownames()

I compared the results and found out that the gene lists were pretty different and shared only 588 of 2000 genes:
shared_variable_genes <- intersect(variable_genes_Seurat, variable_genes_M3Drop)
length(shared_variable_genes)
> 588

My questions are:

Why are the results so different?
Which feature list - shared_variable_genes, variable_genes_Seurat or variable_genes_M3Drop - should I use for dimensionality reduction and clustering?

Thank you!


text:
pip installation within conda env: Error downgrading scikit-learn - WARNING: No metadata
I am trying to get the bioinformatics tool deepbgc up and running. Now I have created a conda environment and have currently installed deepbgc onto it. Now when I run the tool, the program works but it outputs the following warning message as a result:
 /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/sklearn/base.py:306:
 UserWarning: Trying to unpickle estimator RandomForestClassifier from version 0.18.2
 when using version 0.21.3. This might lead to breaking code or invalid results. Use at your own risk.

This led me to believe that the program is trying to utilize an older version of the scikit-learn package (version 0.18.2 instead of the currently installed 0.21.3 version) for this RandomForestClassifer estimator. This pointed me to the next logical solution of downgrading my current version of scikit-learn. So, I ran a simple pip install of the older package version:
pip install scikit-learn==0.18.2

For some reason, the package downgrade failed and outputted the following error:
    error: Command "g++ -pthread -B /EFS/tools/miniconda/envs/deepbgc2/compiler_compat -
Wl,--sysroot=/ -Wsign-compare -DNDEBUG -g -fwrapv -O3 -Wall -fPIC -
I/EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/numpy/core/include -
I/EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/numpy/core/include -
I/EFS/tools/miniconda/envs/deepbgc2/include/python3.7m -c
sklearn/cluster/_dbscan_inner.cpp -o build/temp.linux-x86_64-cpython-
37/sklearn/cluster/_dbscan_inner.o -MMD -MF build/temp.linux-x86_64-cpython-
37/sklearn/cluster/_dbscan_inner.o.d" failed with exit status 1

      [end of output]

  note: This error originates from a subprocess, and is likely not a problem with pip.
  WARNING: No metadata found in /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages
  Rolling back uninstall of scikit-learn
  Moving to /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/scikit_learn-0.21.3.dist-info/
   from /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/~cikit_learn-0.21.3.dist-info
  Moving to /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/sklearn/
   from /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/~klearn
error: legacy-install-failure

× Encountered error while trying to install package.
╰─> scikit-learn

note: This is an issue with the package mentioned above, not `pip`.

Currently, I have Python version 3.7 installed on the environment and version 22.1.2 of pip. Any guidance appreciated, thanks!

pip and conda don't mix and good advice is don't ever mix them. Thus to avoid this delete your old environment, create a new conda environment, or use a conda environment where pip has not been used ...
conda install -c anaconda scikit-learn=0.18.2

If you are using pip it is recommended to use it only with virtualenv venv. If there is a dependency issue delete the environment and start over is much easier.
The rationale is for not mixing them is:

conda installs lots of stuff not just Python and it gets the env (environment) right for this purpose
pip only installs Python within a Python environment

The simple answer is
conda install -c bioconda deepbgc

conda will then install the correct version of scikit-learn automatically. If doing this it important it is in its own (new) conda environment however.
If the above doesn't work further stuff can be done, but the 'simple answer' should work.

@Joe_G says its hanging ... on solving environment just leave it ... you will then see why,
The following NEW packages will be INSTALLED:

  _tflow_select      pkgs/main/osx-64::_tflow_select-2.3.0-mkl
  absl-py            pkgs/main/noarch::absl-py-0.15.0-pyhd3eb1b0_0
  appdirs            pkgs/main/noarch::appdirs-1.4.4-pyhd3eb1b0_0
  astor              pkgs/main/osx-64::astor-0.8.1-py37hecd8cb5_0
  biopython          pkgs/main/osx-64::biopython-1.76-py37h1de35cc_0
  blas               pkgs/main/osx-64::blas-1.0-mkl
  bottleneck         pkgs/main/osx-64::bottleneck-1.3.5-py37h67323c0_0
  brotli             pkgs/main/osx-64::brotli-1.0.9-hca72f7f_7
  brotli-bin         pkgs/main/osx-64::brotli-bin-1.0.9-hca72f7f_7
  c-ares             pkgs/main/osx-64::c-ares-1.18.1-hca72f7f_0
  ca-certificates    pkgs/main/osx-64::ca-certificates-2022.07.19-hecd8cb5_0
  certifi            pkgs/main/osx-64::certifi-2022.9.14-py37hecd8cb5_0
  cycler             pkgs/main/noarch::cycler-0.11.0-pyhd3eb1b0_0
  dataclasses        pkgs/main/noarch::dataclasses-0.8-pyh6d0b6a4_7
  deepbgc            bioconda/noarch::deepbgc-0.1.23-py_0
  fftw               pkgs/main/osx-64::fftw-3.3.9-h9ed2024_1
  fonttools          pkgs/main/noarch::fonttools-4.25.0-pyhd3eb1b0_0
  freetype           pkgs/main/osx-64::freetype-2.11.0-hd8bbffd_0
  gast               pkgs/main/noarch::gast-0.5.3-pyhd3eb1b0_0
  giflib             pkgs/main/osx-64::giflib-5.2.1-haf1e3a3_0
  grpcio             pkgs/main/osx-64::grpcio-1.42.0-py37ha29bfda_0
  h5py               pkgs/main/osx-64::h5py-3.7.0-py37h4a1dd59_0
  hdf5               pkgs/main/osx-64::hdf5-1.10.6-h10fe05b_1
  hmmer              bioconda/osx-64::hmmer-3.3.2-h9722bc1_2
  importlib-metadata pkgs/main/osx-64::importlib-metadata-4.11.3-py37hecd8cb5_0
  intel-openmp       pkgs/main/osx-64::intel-openmp-2021.4.0-hecd8cb5_3538
  joblib             pkgs/main/noarch::joblib-1.1.0-pyhd3eb1b0_0
  jpeg               pkgs/main/osx-64::jpeg-9e-hca72f7f_0
  keras              pkgs/main/osx-64::keras-2.2.4-0
  keras-applications pkgs/main/noarch::keras-applications-1.0.8-py_1
  keras-base         pkgs/main/osx-64::keras-base-2.2.4-py37_0
  keras-preprocessi~ pkgs/main/noarch::keras-preprocessing-1.1.2-pyhd3eb1b0_0
  kiwisolver         pkgs/main/osx-64::kiwisolver-1.4.2-py37he9d5cce_0
  lcms2              pkgs/main/osx-64::lcms2-2.12-hf1fd2bf_0
  lerc               pkgs/main/osx-64::lerc-3.0-he9d5cce_0
  libbrotlicommon    pkgs/main/osx-64::libbrotlicommon-1.0.9-hca72f7f_7
  libbrotlidec       pkgs/main/osx-64::libbrotlidec-1.0.9-hca72f7f_7
  libbrotlienc       pkgs/main/osx-64::libbrotlienc-1.0.9-hca72f7f_7
  libcxx             pkgs/main/osx-64::libcxx-14.0.6-h9765a3e_0
  libdeflate         pkgs/main/osx-64::libdeflate-1.8-h9ed2024_5
  libffi             bioconda/osx-64::libffi-3.2.1-1
  libgfortran        pkgs/main/osx-64::libgfortran-5.0.0-11_2_0_h246ff09_26
  libgfortran5       pkgs/main/osx-64::libgfortran5-11.2.0-h246ff09_26
  libpng             pkgs/main/osx-64::libpng-1.6.37-ha441bb4_0
  libprotobuf        pkgs/main/osx-64::libprotobuf-3.20.1-h8346a28_0
  libtiff            pkgs/main/osx-64::libtiff-4.4.0-h2ef1027_0
  libwebp            pkgs/main/osx-64::libwebp-1.2.2-h56c3ce4_0
  libwebp-base       pkgs/main/osx-64::libwebp-base-1.2.2-hca72f7f_0
  llvm-openmp        pkgs/main/osx-64::llvm-openmp-14.0.6-h0dcd299_0
  lz4-c              pkgs/main/osx-64::lz4-c-1.9.3-h23ab428_1
  markdown           pkgs/main/osx-64::markdown-3.3.4-py37hecd8cb5_0
  matplotlib-base    pkgs/main/osx-64::matplotlib-base-3.5.2-py37hfb0c5b7_0
  mkl                pkgs/main/osx-64::mkl-2021.4.0-hecd8cb5_637
  mkl-service        pkgs/main/osx-64::mkl-service-2.4.0-py37h9ed2024_0
  mkl_fft            pkgs/main/osx-64::mkl_fft-1.3.1-py37h4ab4a9b_0
  mkl_random         pkgs/main/osx-64::mkl_random-1.2.2-py37hb2f4e1b_0
  munkres            bioconda/noarch::munkres-1.0.7-py_1
  ncurses            pkgs/main/osx-64::ncurses-6.3-hca72f7f_3
  numexpr            pkgs/main/osx-64::numexpr-2.8.3-py37h2e5f0a9_0
  numpy              pkgs/main/osx-64::numpy-1.21.5-py37h2e5f0a9_3
  numpy-base         pkgs/main/osx-64::numpy-base-1.21.5-py37h3b1a694_3
  openssl            pkgs/main/osx-64::openssl-1.1.1q-hca72f7f_0
  packaging          pkgs/main/noarch::packaging-21.3-pyhd3eb1b0_0
  pandas             pkgs/main/osx-64::pandas-1.3.5-py37h743cdd8_0
  pillow             pkgs/main/osx-64::pillow-9.2.0-py37hde71d04_1
  pip                pkgs/main/osx-64::pip-22.2.2-py37hecd8cb5_0
  prodigal           bioconda/osx-64::prodigal-2.6.3-ha5712d3_4
  protobuf           pkgs/main/osx-64::protobuf-3.20.1-py37he9d5cce_0
  pyparsing          pkgs/main/osx-64::pyparsing-3.0.9-py37hecd8cb5_0
  python             pkgs/main/osx-64::python-3.7.7-hc70fcce_0_cpython
  python-dateutil    pkgs/main/noarch::python-dateutil-2.8.2-pyhd3eb1b0_0
  pytz               pkgs/main/osx-64::pytz-2022.1-py37hecd8cb5_0
  pyyaml             pkgs/main/osx-64::pyyaml-6.0-py37hca72f7f_1
  readline           pkgs/main/osx-64::readline-8.1.2-hca72f7f_1
  scikit-learn       pkgs/main/osx-64::scikit-learn-1.0.2-py37hae1ba45_1
  scipy              pkgs/main/osx-64::scipy-1.7.3-py37h214d14d_2
  setuptools         pkgs/main/osx-64::setuptools-63.4.1-py37hecd8cb5_0
  six                pkgs/main/noarch::six-1.16.0-pyhd3eb1b0_1
  sqlite             pkgs/main/osx-64::sqlite-3.39.3-h707629a_0
  tensorboard        pkgs/main/osx-64::tensorboard-1.14.0-py37h80053f4_0
  tensorflow         pkgs/main/osx-64::tensorflow-1.14.0-mkl_py37h085be34_0
  tensorflow-base    pkgs/main/osx-64::tensorflow-base-1.14.0-mkl_py37h5a24fda_0
  tensorflow-estima~ pkgs/main/noarch::tensorflow-estimator-1.14.0-py_0
  termcolor          pkgs/main/osx-64::termcolor-1.1.0-py37hecd8cb5_1
  threadpoolctl      pkgs/main/noarch::threadpoolctl-2.2.0-pyh0d69192_0
  tk                 pkgs/main/osx-64::tk-8.6.12-h5d9f67b_0
  typing_extensions  pkgs/main/osx-64::typing_extensions-4.3.0-py37hecd8cb5_0
  werkzeug           pkgs/main/noarch::werkzeug-2.0.3-pyhd3eb1b0_0
  wheel              pkgs/main/noarch::wheel-0.37.1-pyhd3eb1b0_0
  wrapt              pkgs/main/osx-64::wrapt-1.14.1-py37hca72f7f_0
  xz                 pkgs/main/osx-64::xz-5.2.6-hca72f7f_0
  yaml               pkgs/main/osx-64::yaml-0.2.5-haf1e3a3_0
  zipp               pkgs/main/osx-64::zipp-3.8.0-py37hecd8cb5_0
  zlib               pkgs/main/osx-64::zlib-1.2.12-h4dc903c_3
  zstd               pkgs/main/osx-64::zstd-1.5.2-hcb37349_0

Proceed ([y]/n)?

Its an enormous data pipeline. Everything is there, thats why it takes so long.

Hi, thanks for your response. I went ahead and created a new conda environment and tried to run the deepbgc install on there with conda instead of pip. The problem now is the install is getting stuck at the "Solving Environment" phase. For whatever reason, It doesn't seem to make it past that point.

Hi @Joe_G I've updated the post, so did it work?

I've found mamba to be a life-saver for getting around conda's slowness with large sets of packages.  Its mamba command can act as a drop-in replacement for conda for most purposes and most importantly it solves dependencies much faster.

Unfortunately it never got past the solving environment stage. I was installing it on a remote server and the process took so long that the my user was disconnected from the server before it could finish. I'm going to try again thinking maybe it was an issue with the server. @Jesse I may try Mamba that sounds like it could be helpful. Is that just a package installation? As mentioned before, the warning message wasn't code breaking so I may just be living with it for the time being. I do appreciate the guidance, thanks!

Thanks for your feedback @Joe_G . Personally conda has been the most stable environment and package management system I've used. Upvotes are always appreciated BTW.

Will do @M_ Thanks again!

Yes mamba itself is just an anaconda package so it's a simple install.  You have to install it into the base environment though (something like conda install mamba -n base -c conda-forge to get the latest one from conda-forge.)


text:
Bulk protein annotation with Prokka
Can prokka be used to annotate a set of protein sequences (so bulk protein annotation), or does it require an entire genome for the annotation?
I've been looking through the documentation but haven't found a clear answer to this. I tested it, and the annotation of my protein sequences worked, but I wanted to check if the program expected a genome and that I can trust the output.
Thanks!

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

Prokka doesn't need whole genomes, I've definitely used it for smaller contigs before so I think you're fine!

Thanks for the reply, that's great!


text:
How to get started with louvain/leiden algorithm with UMAP in R
A collegue of mine recently suggested to try the louvain algorithm for clustering multiplex cytometry data. However, implementations of louvain are kind of rare in R. To my knowledge the only stand-alone implementation is included in the igraph package.
So far i tried to use the adjacency matrix that the umap function returns (uwot package) to create a graph using igraphs graph_from_adjacency_matrix function. I can run the louvain algorithm on the graph, but the result is always a few thousand clusters with a hand-full if cells. changing the resolution parameter does not change anything.
If i remember correctly, Seurats findClusters function uses louvain, however i don't want to use PCA reduction before clustering, which is requiered in Seurat to find clusters.
Does anybody know of a different implementation of louvain or leiden algorithm that or an easier way to use igraphsfunctions? I obviously don't have to use the matrix umap returned, but what would be a better approach to make my data "louvain-able"?

Both conda and PyPI have leiden clustering in Python which operates via iGraph
conda install -c conda-forge leidenalg

pip install leiden-clustering

Used via
import leidenalg as la
import igraph as ig

Example output

The docs are here. You will not need much Python to use it.
There is an entire Leiden package in R-cran here
Python implements Rs Leiden package here,
conda install -c conda-forge r-leiden


text:
Why are genetic algorithm chromosomes represented using 8 bit binary?
I am learning about genetic algorithms.
Why chromosomes are represented using 8 bit binary.
This arose when I apply crossover techniques to the bits.
Please help me to find the answer.Thanks in advance !

Hi @Julie I do understand you originally posted on Stackoverflow and were directed here. This is a good question, however it is a statistics problem because the chromosomes are not a physical chromosome but a name given to the algorithm design. https://stats.stackexchange.com is better suited to this question.

I agree that it's not biology, but I consider genetic algorithms to be a programming problem, not a stats problem (and have been asked and answered on SO previously, e.g.). Seems to me like StackOverflow would have been an appropriate place to ask this question.

Not always. The representation of chromosomes in genetic algorithms is often dependent on the specific problem to be solved and number of variables to be found through the optimization process.
Not only you can chose a 16-bit representation, but integer-based representations are also used as mentioned here.
At the computation level, using a binary representation usually facilitates cross over operations on the bit arrays. 8 bits are compact while already give you a range of 0-255 different values for a given variable, thus provides a first option to go unless you need a larger value space.


text:
reliable and popular HLA imputation tools
Could someone recommend me popular and reliable tools for HLA imputation?
I'm starting my master thesis project and I will infer the HLA region of 2000 sample in a case-control study; I searched for benchmarks on PubMed but I really wanted to now if someone experienced in this kind of analysis can help me
ty

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.

You at the very least have to specify what kind of data you are working with.

This sounds like a question that is opinion based, without a clear best answer, and an answer that is likely to change in the future. Such questions are difficult to handle  with the stackexchange format. May I suggest you try https://www.biostars.org/ ?


text:
How can I get the HGVS notation of a mutation from a FASTA?
In an old article, they gave a table of mutations found in a gene and its positions which points to another article containing the FASTA sequence.
This FASTA is the precursor sequence of the gene. How can I know which transcripts those mutations belong to? And how can I obtain their HGVS format?
For example, C-7111-G should have an HGVS notation like c.XXXC>G.

What's the gene? There are tools you can try, but we'd need to know the gene first. And no matter what you do, this will be a manual process since the old FASTA sequence you have is likely inaccurate anyway.

The gene described in the article is POMC. What I did is the following: first, I copied the same fasta sequence from the old article, then I pinpointed the exact location of the mutations from the table. Second, I copied a couple of unique nucleotides from the left and right of the mutations. Third, I matched those sequences from the actual gene fasta sequence to locate the real nucleotide, and I got its position in the current hg38 build. Then the rest is simple, I searched in the databases what mutations are found in those emplacements.

Yeah, that's probably the best thing you can do. I was going to shamelessly plug the product of the company I work for, https://varsome.com, which is free and can help you translate variant notations, but if all you have is the position relative to an old fasta sequence and don't know the actual transcript, we can't really help much more. You'll be able to get the proper HGVS notation easily though. Same for https://variantvalidator.org/ with which I am not affiliated and which is a great tool for validating HGVS.

Thank you @terdon. I will check those websites. Cheers!