text:
Extracting positions from PAF files in order to extract sequences from a Fastq file with Python
I have used Minimap2 to create a paf file by aligning a Fastq file against itself. Now from this Paf file I can see where the reads overlap, and I want to take these positions, and use them to extract the entire sequences between these positions from the paf file so like this:
If read A starts on position 72 and ends on position 2041 and overlaps here with read B which starts on position 190 and ends on position 2741, I would like to extract these from the Fastq and merge into a contig. But I don't really know how to do that.
query_name,query_length,query_start,query_end,strand,target_name,target_length,target_start,target_end,residue_matches,alignment_block_length,mapping_quality,tp,cm,s1,dv,rl
read=73,2677,81,2674,+,read=86140,3729,190,2779,2082,2627,0,S,4,2,0,6
read=73,2677,89,2674,-,read=78530,52349,15309,17890,1977,2623,0,S,4,1,0,6
read=73,2677,95,2674,-,read=8246,38338,32497,35071,1976,2623,0,S,4,1,0,6
read=73,2677,95,2674,-,read=82661,13059,8472,11069,1960,2636,0,S,4,1,0,6
read=73,2677,100,2669,+,read=79559,11152,5925,8498,1928,2617,0,S,4,1,0,6
read=73,2677,89,2674,+,read=75579,16131,10030,12576,1924,2609,0,S,4,1,0,6
read=73,2677,100,2669,+,read=92760,9659,7065,9617,1904,2613,0,S,4,1,0,6
read=73,2677,89,2674,+,read=95665,81160,4544,7111,1866,2624,0,S,4,1,0,6
read=73,2677,81,2660,+,read=75036,31663,27966,30529,1841,2619,0,S,4,1,0,6

This is a snippet from the PAF file. It's comma delimitted, and below are some of the reads from the Fastq file, the reads vary a lot in length.
@9a65dba6-05c9-49d7-ac0e-2ede30be6036 runid=7f8eec05f1beed68c29c091f78bf84807beb3275 read=39 ch=1736 start_time=2021-01-27T22:53:24Z flow_cell_id=PAG49596 protocol_group_id=210105Chi sample_id=D21-130 barcode=barcode01
AGTATACTTCGTTTCAGTTACGTATTGCTAAGGTTAACACAAAGACACCGACAACTTTCTTCAGCACCTGGTGCTGAAGAAAGTTATCGGTATTTCTTTGTGATCAGCCGTAG
+
%))(&&''&((AA8>33??7++)'*+$$55(1@BFFF>?>=>684.-8:2)&.,11168=>,*/13;;;348:A>;7032)2312%*,&')(001123)),+&,*$%#$%$$$
@edaf3d67-a971-43ed-be5b-3816b985d53f runid=7f8eec05f1beed68c29c091f78bf84807beb3275 read=136 ch=1940 start_time=2021-01-27T22:54:06Z flow_cell_id=PAG49596 protocol_group_id=210105Chi sample_id=D21-130 barcode=barcode01
ATTACTTCGTTCAGTTACGTATTGCTAAGGTTAACACAAAGACACCACGACAACTTTCTTCGCACTATGGTGTTCGCTATTATCAGATTCCTGATGAAGAAGAACTAATGAAATTGGTTTCTATTATTTATTGACTTCGTTCTATAATAGAGCAAAGTTTTCAGCTGGAGAAAATTGGTGAAGAAAAAGCTGAGATTACTGCAAAAGCAACTTCTCTCGTGGTAAATCTGTCACCTCTGCAGAGCATTATCTGAAAGAATGAGGAACTTTCCTGAAGTTATTGTTCAATATGGTTGCTCTTCAAAATTGCTAAGCCAGAACAGCAAGATTGATAACATTCATGCATTAGAAGACGCTCTTTATAGTAAGTTGGTTGGCGTGGTACAGTAGACTGTATTAGTGAATATGATGGAGAATTAGCAGTTATTGACTTCAAGACTTCTAAAGCACCTAAGACCGAAAATGGATTGATGGTTATTTTGCAGCAATGGCATGTATTTATTCATTAGTTCATTGGCATGTGCTAAGAAACTTGTTGTATAATAATGGCATGTGAAGATGGAACCTGTAAAGTTTATGAAGGCGTAATAGAAACATCAACATCAGCAAATTACTCGTACATTACATTAGAAACTTTTTGGATTTTCAACTAGAATTAAATGGAAAATGAATTAACCAAAGCACTAGATAAAAAATTCATGAATGCTGCCAAAAAGTTTTCTCTTGAAGTAGAGAAAATTGTATTGGAAGAAAAAAGTGAATTATATTGATGCTATTGTTCTATTTTGTGGAAAAATAGTATTGAAGTGGACTCTGTTACTAAGTTGATTTCTAAACCTTTAAAGGAAAAATTAAAGGCAGATACTGTAATTTGAACTTCATGAAGAACCTCGTGCCGACTACCCATATAAATACTAGTATCTAGGAGGAATTATGTCTGATTTTTTGAATCTGAACAAGTTCCAAGAATCTATAAAGACATCAACAGATGCAAGAAAGAGATATATTCTAAAGTCTTGTTTGCTTTTTGACATTGCCTCAAAGAAAATGGATCATGTTGATTTATTAGAACAACTTCTAGAGAAGCAAAAGAAACTTCTACATGAGATTGGCTTCTGACGATCCTCGTGCTAAAGATATGATGGAGCAAATTCATCAATCTGCTCAGTTGATTGGATTCCCCTAAGGATTCATCACCTGATGTTTTATTTGAGAATATGAACCTTGGTCTTATTAGATTTGCTTGACATGGGAGCACTGTCGCCCTATAATAGACCAGTACACAAAAGCCCAAATCTAATTACAGAAATTTGGAAATCTATGTCTTTTGCATCACTAAAAAAACAGTCCTCTCTTGGATCTCTAACCGCAAAGTTGGTTAAAGAAGTAGAAGACCAATAAAGCAGGTAACGGGCTGATGAGCGTCTTTGACCAGAGAGTTGATAAAAAGCATGGAGCAGTCTGTTGCAGTTATTCGTTTCTTCCAGCACCAGATAAGAAGGATTCTTCCTTGGGCAAAACTTTATTCTCATGCATTTCAAGGATGAGTGGATGGTACGTGAAAATAACTACTGTTGGTAGTAAGGATCCTGTTTCAGAGTACAATACGGTCTCATGGATTCTGGTTCAGATGCCAATAAGCAGATCGCTCCTTCAATCAAAGGCGTAAGCTTTTCTTACTATAGCAACATCTATGTTGTAGAGGATATACAAATCCGAATGGAAGGAAAAATATGTATTTACAAATTTGGTAAGAAGATCTTGACAGATCATCGCCAATAATACAACCTGAGTTCAGGGATGAGTCACCAATCAATCCATTGATTTCACAA
+
%%&-79888F?9.-90/5400/0(&24105EEH?:39<7++999:8'+23695,,,..77*%&'*.'%).&)882))2*220.7,3CIHA7476935-2734D;:;8+--1;003))%'&$%$$&#%&&&&(-57:;44>:;B=<A4//.%+452,-4('(#$##%&++&,.9;???:CEHJ9GEA@</10.:9<HH2<<;>IE10,45*79:1(3%1,3-*)'''*'''31*++40*$$$$)'689<?A5668:)4<=220#''%$$$$$%%%.++2/%#%$'(#+*)(&#&')#$')*+'>><::<65<:,-:;8)&((--4;9>=ED?<<==;;886)6689<;&&)&%$"$-55;>=+9:29;)'('')#%%%$$&&&%**+112&).739:@=-0.-.290<>>@DABCFGFGABBDEIHD>>02'''(/-6977>64;62,+/33/&(#**+,.7?:--07;:74'')4.''((&%$((*100.)%%$##$##$#$%#%#%$$%*$&#$$$$#$%'&%((4-;<79992..++2589=B?;;<;=?@>@DE=FAE>H@?===ADBA=HK@>=;<+(830-1--($*+%$)%&$'#(&%%$$%.5569CH@=>==A@A=CGHJKKNJFEKGB0'(1/HLJA<<'4:7/?>99ABE??@CJNG==>3?EJFCCROOFEGD>EHJLWTQMMJIMNJ?@CFBBC@>=>9:2/.1>CEF?GAFDC;@5::&8(BBIFJHBCEFF;571*II?34=6;=BFFG:55=;??7.+*,56@KLFA@@@C<,(+/0,>E734867526/;<DB<DEDDG?=F:>D@?CI@EB@CMHEF@<@B@A>80,444=AC7781)%%&&%(%$#$&#&''(6-.35''7'%(()&)&()*>74<8<.)038./0.&@BDICA=>KE@<<CC><7<;8747?=>;66<:;9ANFED4%4762/0//3;<8:DB8:=8+;<>:E?@D@/59<<956/13,=88<B6@?=',:8@><89>@;)()+%.*&''())&./0/,%'(&&''(2263'''%07A<;<<C?:?===B<8:JD6;===?A@CB//2/.(%$14333@<1.)31,5643+.:;89=-/9:+&)-./022149892163/+;56>=;,'''2360-..;ABC@@98;<@>E=<?>AA83C>==C?*.*FGE@;;3-.---AD9;<.;89=??;;?HA/$$'81*.3+(%#$%($$#%%$'&$&$%&%*%*,92'')7;8:9;?;9=???8:@76;>=A>>A<;>8=9>;==A9<8?==;(1+).1AC=ECDFF501&*)&#(5258EAB>@&/285>@?3.;@@<:=<A@><;;3009;5+556<AA?78=7C;>?;67+7:;75>;,/92./,'++.//(()/59<@>>CA8,&&8:894%%(.%/01//.0.,+',/'&,/10)%$$(/,&&$.,***#%%%&&$##$$%%##$'&&')26&1331:;A:CD?A=3666:/.300(,*/#'#'+#1;6444-+(+35<;89<7257,@@>1.++DDA4*&(*'*)*7<<<=873)#%%$#$--(''%(*+7968>6.41-7888:>B:;==FDC:6998;;8:+'&$$'+))'#(%+%')$(+2<?<56:3116708;8%445433.474/-1+--5465'&&%$&++,,088:5:;95;/&*::@99=A::;45;?)$$&+,,&$$#$%00--.0(),''(&''''&&&$$%##$%%&$'224(12;.001+$*)(()'%&%##$$&'%$%$#"$))$0+)'8;79:;93=A?;&%136-,,/;.+*3;<>8?>=:<9<<58?98.#$&%

Another thing that is currently confusing me about my Paf file is the target length and query length which seems to make no sense. It's not the length of the overlap nor the length of the entire segment or read. Can anybody shed some light here on the Paf file and how to extract these positions in Python? Please note that I do not want to use external tools for this, only Python and Python libraries.

The default output from minimap2 is tab-separated, not comma-separated, and doesn't have headers. Could you please show how you produced that CSV file (either code or a command line would be helpful)?

Your example CSV output has a query name of "read=73" for all lines,  but the fastq file doesn't have any sequences with that format for read names. It has channel read sequences, but those are different (and also don't match). Could you please add that read into the fastq example (either by using the read uuids in the CSV file, or renaming the sequence IDs in the example fastq file to match the csv file)?

@gringer the read names do have read=NN as the 3rd field, perhaps that is being picked up by the tool? Or are those the "channel reads" you mentioned?

Each flow cell has a fixed number of sequencing channels; reads from each channel are numbered sequentially, and this is the read=NN metadata you see. This is not unique: it's possible to have the same read number from different channels. The channel read sequence should not be used as the read ID; it's better to use the UUID for that (i.e. follow the standard process for read names and use the value after @).

The CSV file was made by an external script, the headers were added that way, so i could understand what was going on, the internal structure of the PAF is still the same.

Please show the command line or script that did the transformation from the PAF file; the read and reference names have been changed from what I would expect from a direct transformation from a PAF file.


text:
sra toolkit fastq-dump comand not terminating
I ran the following command trying to download a fasta file from NIH but the command doesn't seem to terminate and I get no output file in my repository:
fastq-dump SRR10134980

I have installed SRA Toolkit following the instructions on this page.

Do you get the same behavior when you use fasterq-dump? Note that fasterq-dump is the successor to the older fastq-dump tool.

That's an 11GB file and fastq-dump is not a speedster. Without more info the likely answer is that it keeps downloading which might take several hours on  poor connection. You can get fastq download links from sra-explorer.info -- That is more stable than this tool.

source: Stack Exchange

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

sort: votes


text:
Nextflow Pipeline: Unexpected Input "{"
I'm attempting to run a nextflow pipeline and have begun getting unexpected output for a few of the processes in my pipeline. I believe this to be a syntax error, but I am not sure where to correct it. Any help is finding where my error is located is appreciated.

#!/usr/bin/env nextflow
 nextflow.enable.dsl=2
 import groovy.io.*

  process hisat2Alignment {
    fastq1 = NS.1389.002.NEBNext_dual_i7_*_R1.fastq.gz
    fastq2 = NS.1389.002.NEBNext_dual_i7_*_R2.fastq.gz
    input:
    tuple path("fastq1"), path("fastq2")

    output:
    tuple path("*.sam"), val(fastq1.baseName.substring(0, fastq1.baseName.indexOf("_R1")))
    script:
    """
    hisat2 -p 16 -x genome -1 ${fastq1} -2 ${fastq2} -S mapped_${baseName}.sam
    """

 }

process Stringtier {
    publishDir "${params.baseDir}/assembled/$baseName", mode: "move"
    env {
     PATH = $PATH:/data/selena/stringtie-2.2.1
    }

    input:
    tuple path(mapped_bam), val(baseName)
    output:
    tuple path("*.gtf"), val(baseName)

    script:
    """
    stringtie -e -G ${params.annotation} -o ${baseName}.gtf ${mapped_bam}

    """
}

 workflow {
    samples = new File(params.configFilePath).readLines()*.split('\t');
    samples  = Channel.fromList(samples);
    hisat2Alignment(samples);
    assembled = Stringtier(samples)

}

Nextflow does not have an env process directive. If stringtie is not already in your $PATH, you'll need to append it before running nextflow, for example:
export PATH="${PATH}:/data/selena/stringtie-2.2.1"
nextflow run de-pipeline.nf

However, you may prefer to avoid local tool installations altogether, and use either the conda or container process directives. The latter will however require a container runtime, like Docker or Singularity.
You may also consider including a 'sample name' or some sort of sample identifier in your config file, which I suspect is just a simple tab-separated file of filenames. This would greatly simplify your workflow. For example:
params.samplesheet = './samples.tsv'
params.genome = './genomes/GRCh38/GRCh38'
params.gtf = './Homo_sapiens.GRCh38.gtf'

genome_name = file( params.genome ).name
genome_path = file( params.genome ).parent

params.outdir = './results'

process HISAT2 {

    tag { sample }

    publishDir "${params.outdir}/HISAT2/${sample}", mode: 'copy'

    cpus 16

    input:
    tuple val(sample), path(reads)
    path genome_dir

    output:
    tuple val(sample), path("${sample}.sam")

    script:
    def (r1, r2) = reads

    """
    hisat2 \\
        -p ${task.cpus} \\
        -x "${genome_dir}/${genome_name}" \\
        -1 "${r1}" \\
        -2 "${r2}" \\
        -S "${sample}.sam"
    """
}

process SAMtools_sort {

    tag { sample }

    input:
    tuple val(sample), path(bam)

    output:
    tuple val(sample), path("${sample}.sorted.bam")

    """
    samtools sort \\
        -o "${sample}.sorted.bam" \\
        "${bam}"
    """
}

process StringTie {

    tag { sample }

    publishDir "${params.outdir}/StringTie/${sample}", mode: 'copy'

    input:
    tuple val(sample), path(bam)
    path gtf

    output:
    tuple val(sample), path("${sample}.gtf")

    """
    stringtie \\
        -e \\
        -G "${gtf}" \\
        -o "${sample}.gtf" \\
        "${bam}"
    """
}

workflow {

    gtf = file( params.gtf )

    Channel
        .fromPath( params.samplesheet )
        .splitCsv( sep: '\t' )
        .map { sample, fq1, fq2 ->
            def reads = tuple( file(fq1), file(fq2) )

            tuple( sample, reads )
        }
        .set { samples }

    HISAT2( samples, genome_path )

    SAMtools_sort( HISAT2.out )

    StringTie( SAMtools_sort.out, gtf )
}

Where the contents of the input sample sheet (e.g. samples.tsv) looks like:
sample_A    /path/to/reads/A_R1.fastq.gz    /path/to/reads/A_R2.fastq.gz
sample_B    /path/to/reads/B_R1.fastq.gz    /path/to/reads/B_R2.fastq.gz
sample_C    /path/to/reads/C_R1.fastq.gz    /path/to/reads/C_R2.fastq.gz

Results:
$ nextflow run main.nf
N E X T F L O W  ~  version 22.04.0
Launching `main.nf` [adoring_poisson] DSL2 - revision: 43e8c59949
executor >  local (9)
[60/1cba2f] process > HISAT2 (sample_A)        [100%] 3 of 3 ✔
[0a/931a07] process > SAMtools_sort (sample_A) [100%] 3 of 3 ✔
[6d/7689a2] process > StringTie (sample_A)     [100%] 3 of 3 ✔

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


text:
Simes method in R
I have data that looks like this: 3 column SNPs their gene based on Annovar and a p-value for every SNP. What I would like is to aggregate the p values for every gene.
snps <- data.frame(
  snp_id = c("rs1", "rs2", "rs3", "rs4", "rs5", "rs6", "rs7", "rs8"),
  gene = c("gene1", "gene1", "gene1", "gene1", "gene2", "gene2", "gene2", "gene2"),
  pvalue = c(0.7703884, 0.9648540, 0.9648540, 0.9648540, 0.54, 0.03, 0.03, 0.8)
)

what I need to do is :

For each gene, sort the p-values in ascending order.
For each p-value in the sorted list, multiply it by the number of SNPs in that gene.
Divide the result from step 2 by the rank of the p-value in the sorted list (i.e., its index in the list).
Take the minimum value from step 3 as the new p-value for that gene.
Adjust the p-values using the false discovery rate (FDR) correction method.

This is the formula for each gene :combined p-value $\min(n p_{(i)}/i)$ where $p_{(1)}, …, p_{(n)}$ are the ordered p-values.
I tried this code  with the usage of simes.test package in R:
library(mppa)

x=c(0.77,0.964,0.964,0.964)
simes.test(x)

The answer is 1.285. This is the source code of the simes.test function:
function (x, returnstat = FALSE)
{
    r = rank(x)
    T = min(length(x) * x/r)
    if (returnstat)
        c(T, T)
    else T
}

How is it possible for the pvalue to be larger then 1? What is the problem? I would be happy if you could share your knowledge in this situation.

What is a p-value for an SNP?

It can happen it seems. detailed discussion of this topic can be found here http://stat.wharton.upenn.edu/~steele/Courses/956/ResourceDetails/MultipleComparision/Writght92.pdf


text:
Convert Abricate output tsv file to gff3 format
I have a tsv file generated from abricate (https://github.com/tseemann/abricate). I need to convert them to gff3 format with certain columns retained, certain columns reordered, while other columns deleted.
We are trying to use these gff3 files for downstream applications and for piping into other applications. However, we could not solve it.
Below are examples of my tsv files, what possibly we may need to do, and desired output files in gff3 format.
Any help will be much appreciated.
Input tsv file:
FILE SEQUENCE START END GENE COVERAGE COVERAGE_MAP GAPS %COVERAGE %IDENTITY DATABASE ACCESSION PRODUCT
UBird_Cyou_D3.fna BJCZ01000001.1 1866608 1867417 cdtB 1-810/810 =============== 0/0 100 90 vfdb CAD48850 (cdtB) cytolethal distending toxin B [CDT (VF0185)] [Escherichia coli O157:H str. 493/89] UBird_Cyou_D3.fna BJCZ01000001.1 1867414 1868190 cdtA 1-777/777 =============== 0/0 100 90.61 vfdb CAD48849 (cdtA) cytolethal distending toxin A [CDT (VF0185)] [Escherichia coli O157:H str. 493/89] UBird_Cyou_D3.fna BJCZ01000001.1 2245186 2246238 ompA 1-1041/1041 ========/====== 1/12 100 94.11 vfdb AAF37887 (ompA) outer membrane protein A [OmpA (VF0236)] [Escherichia coli O18:K1:H7 str. RS218]
What we may need to do (there may be other ways too, I am not sure):
Row 1 (always starts with "#") - Need to replace with the string "##gff-version 3"
Col 1 - get rid of ".fna" and retain other data
Insert new Col - print the string from Col 11 for all rows
Col 2 - get rid of entire column
Insert new Col - print "CDS" for all rows
Col 3 - retain data
Col 4 - retain data
Insert new Col and print "." for all rows
Insert new Col and print "+" for all rows
Insert new Col and print "0" for all rows
Col 5 to Col 10 - get rid of all these columns and data
Col 11 - delete column
Col 13 - retain data except "(", ")", "[", "]"
Add new Col - Starting with "ID=" followed by the string taken from Col 1 and a underscore added (for the example data "UBird_Cyou_D3_") alongwith numerals starting from 1 and incrementing by 1. This column data needs to be appended by "product=" followed by data from the corresponding row of the modified Col 13. The separator between ID string and product string should be ";". After completion this column should be like "ID=UBird_Cyou_D3_1;product=cdtB cytolethal distending toxin B CDT VF0185 Escherichia coli O157:H str. 493/89"

Desired final output (*.gff3) considering the example data:
##gff-version 3
UBird_Cyou_D3 vfdb CDS 187 756 . + 0 ID=UBird_Cyou_D3_1;product=cdtB cytolethal distending toxin B CDT VF0185 Escherichia coli O157:H str. 493/89

Hi and welcome to the site! Please [edit] your post and ensure only the actual data are formated as code. We need this in order to be sure we have your data correctly. The rest, your proposed algorithm should be in normal formatting, but your input file and desired gff need to be code formatted.

Is your expected output incomplete? Do you expect four lines of output (including the header)? Are your START and END positions correct?

The awk code is working. However, we need to add little more to take this to maturity. Posted here (https://bioinformatics.stackexchange.com/questions/20761/improving-conversion-of-abricate-tsv-file-to-gff3-file). Any help?

Here's one way using awk, that I think fulfills the requirements. It adds each of the column names (on the first line) to an array to make accessing each of the fields a bit easier. This approach isn't strictly necessary, but it does make for a more readable solution in my opinion. With the following in a file called tsv2gff3.awk:
BEGIN {
    FS=OFS="\t"
    c=1
}

NR==1 && sub(/^#/, "") {
    for(i=1; i<=NF; i++) {
        a[$i]=i
    }

    print "##gff-version 3"
    next
}

{
    sub(/\.fna$/, "", $a["FILE"])
    gsub(/[\(\)\[\]]/, "", $a["PRODUCT"])

    id = "ID=" $a["FILE"] "_" c++
    product = "product=" $a["PRODUCT"]

    print \
        $a["FILE"],
        $a["DATABASE"],
        "CDS",
        $a["START"],
        $a["END"],
        ".",
        "+",
        "0",
        id ";" product
}

And the following in file.tsv:
#FILE   SEQUENCE    START   END GENE    COVERAGE    COVERAGE_MAP    GAPS    %COVERAGE   %IDENTITY   DATABASE    ACCESSION   PRODUCT
UBird_Cyou_D3.fna   BJCZ01000001.1  1866608 1867417 cdtB    1-810/810   =============== 0/0 100 90  vfdb    CAD48850    (cdtB) cytolethal distending toxin B [CDT (VF0185)] [Escherichia coli O157:H str. 493/89]
UBird_Cyou_D3.fna   BJCZ01000001.1  1867414 1868190 cdtA    1-777/777   =============== 0/0 100 90.61   vfdb    CAD48849    (cdtA) cytolethal distending toxin A [CDT (VF0185)] [Escherichia coli O157:H str. 493/89]
UBird_Cyou_D3.fna   BJCZ01000001.1  2245186 2246238 ompA    1-1041/1041 ========/====== 1/12    100 94.11   vfdb    AAF37887    (ompA) outer membrane protein A [OmpA (VF0236)] [Escherichia coli O18:K1:H7 str. RS218]

Run using awk -f tsv2gff3.awk file.tsv, results:
##gff-version 3
UBird_Cyou_D3   vfdb    CDS 1866608 1867417 .   +   0   ID=UBird_Cyou_D3_1;product=cdtB cytolethal distending toxin B CDT VF0185 Escherichia coli O157:H str. 493/89
UBird_Cyou_D3   vfdb    CDS 1867414 1868190 .   +   0   ID=UBird_Cyou_D3_2;product=cdtA cytolethal distending toxin A CDT VF0185 Escherichia coli O157:H str. 493/89
UBird_Cyou_D3   vfdb    CDS 2245186 2246238 .   +   0   ID=UBird_Cyou_D3_3;product=ompA outer membrane protein A OmpA VF0236 Escherichia coli O18:K1:H7 str. RS218

Thank you Steve (https://bioinformatics.stackexchange.com/users/10775/steve) this does work.


text:
Why are asymmetric distributions an issue for GSEA?
I'm reading Subramanian et. al's (2005) original GSEA paper. In the paper, the authors mention the following:

We noticed that the use of
weighted steps could cause the distribution of observed ES scores to
be asymmetric in cases where many more genes are correlated with one
of the two phenotypes. We therefore estimate the significance levels
by considering separately the positively and negatively scoring gene
sets

I understand what they're saying, but I simply don't understand why it's a problem. Basically, they're saying that in some cases, most genes in the list we're testing could appear mostly at the top or bottom of the list. But I don't understand why having asymmetric distributions is a problem. What is the issue here?
EDIT: Here are some example GSEA plots I found via a web search with such gene sets. Where is the problem?

One thing I don't get though is why this problem would be introduced through "weighting" of the steps. I feel like it would also be a problem without weighting. The way it is phrased in the paper seems to suggest that it's the weighting that causes an issue.

The problem is that the enrichment score is modified based on whether the sets are positively or negatively associated, which means that an unbalanced set via either associative statistic or gene count (as is almost always the case) will lead to an incorrect total score.
Put in other words, the calculated score from the positive genes alters the calculated score from the negative genes.
Given this, it's a better idea to either separately calculate scores for the positive set and the negative set (this is the easiest approach, and works with existing methods), or (equivalently) recalculate enrichment score modifiers for positive and negative genes separately (i.e. adjust the graph) so that the zero point for the enrichment score crosses the zero point for association.
A common argument against this (when I have mentioned it to others) is that GSEA has always been done on the total gene set, so it makes sense to do it on the total gene set because it has always been done that way. I'm not a huge fan of the status-quo argument.
I alluded to this issue in this previous answer.
With regards to your figures:
a. Zero-point crossover is at an enrichment score of ~0.42, so the negative set will be underemphasised. It's not contributing much anyway (except for the one most negative gene at the end), so wouldn't affect the outcome much.
b. Zero-point crossover is at an enrichment score of ~0.45. Similar situation as in a.
c. Zero-point crossover is at an enrichment score of ~0.32. There aren't really any positive genes with substantial association scores, so no issue.
d. Plot looks to be identical to c, so the same comments apply.
The scores are most affected when they have a lot of area on both sides of the Y axis. The top right plot of this sub-figure is an example of that, showing a sharp-ish tail on the negative side, so the calculated enrichment score for the negative peak is less extreme than it would be if only the negative genes were considered [source]:

The way it is phrased in the paper seems to suggest that it's the weighting that causes an issue.

Yes. If the positive genes are weighted separately from the negative genes, then the bias can be removed. I did briefly attempt a non-parametric weighting for GSEA (based on rank), but got too confused by the maths involved. Someone answered my question about it, but I haven't had to do any GSEA since then, so haven't actually verified if their solution works for me:
https://math.stackexchange.com/questions/3567714/probability-that-the-nth-highest-of-a-set-of-k-numbers-selected-from-m-numbers-h

Ah I believe I see your point now about recalibrating the graph I think. In practice, I'm guessing this would be done by having N - N_H in P_miss be equal to the "missing genes" on the positive side ONLY, for positive genes, and negative side ONLY, for negative genes

Analogously we'd need to modify P_hit so the denominator only has genes on positive side vs negative side


text:
complete the res type in sdf
I use the diffdock (link to diffdock on github) to dock a peptide with protein and this tool give me a ligand in sdf file format. I traslate this ligand sdf format file to a pdb format file and find the pdbfile lack res information. I wonder if tools/code can automaticly complete this file. I need the UNK become the residue types like pdb.This my error pdb file（this structure is a peptide）

HETATM    1  C   UNK     0       9.797  22.501  -6.484  0.00  0.00           C
HETATM    2  C   UNK     0      10.366  21.086  -6.446  0.00  0.00           C
HETATM    3  O   UNK     0       9.550  20.218  -7.135  0.00  0.00           O
HETATM    4  N   UNK     0      11.710  21.084  -6.978  0.00  0.00           N
HETATM    5  C   UNK     0      11.793  20.201  -8.123  0.00  0.00           C
HETATM    6  C   UNK     0      11.151  20.760  -9.346  0.00  0.00           C
HETATM    7  C   UNK     0      11.586  20.031 -10.581  0.00  0.00           C
HETATM    8  C   UNK     0      13.000  20.308 -10.986  0.00  0.00           C
HETATM    9  C   UNK     0      13.397  19.716 -12.309  0.00  0.00           C
HETATM   10  N   UNK     0      12.636  18.523 -12.633  0.00  0.00           N
HETATM   11  C   UNK     0      11.861  18.392 -13.830  0.00  0.00           C
HETATM   12  O   UNK     0      11.837  19.364 -14.611  0.00  0.00           O
HETATM   13  C   UNK     0      11.148  17.131 -14.065  0.00  0.00           C
HETATM   14  C   UNK     0      11.800  16.246 -15.085  0.00  0.00           C

Please don't post images of text. We cannot use them: we can't copy, we can't search. Also, this makes the post heavier for people with slow connections. Instead, copy/paste the text directly into your question and use the formatting tools to format it as code.

open babel add correct res type to pdb when you feed peptide sdf file, use http://www.cheminfo.org/Chemistry/Cheminformatics/FormatConverter/index.html and https://www.litemol.org/Viewer/ to check it out –

An SDF file lacks residue information, including atom name.
Although for the latter, the molFileAlias atomic property is frequently used.
Adding residue info
There are several ways to do this, such as in RDKit (blog post):
from rdkit import Chem
from rdkit.Chem import AllChem

from typing import List

with Chem.SDMolSupplier('.sdf') as mol_store:
     mols: List[Chem.Mol] = list(mol_store)

mol: Chem.Mol
for mol in mols:
    atom: Chem.Atom
    for i, atom in enumerate(mol.GetAtoms()):
        info = AllChem.AtomPDBResidueInfo(atomName=f'ATM{i}',
                                          residueName='LIG',
                                          residueNumber=1,
                                          chainId='X')
        atom.SetMonomerInfo(info)

mol: Chem.Mol
for mol in mols:
    name: str = mol.GetProp('_Name')
    Chem.MolToPDBFile(mol, f'{name}.pdb')

Notes:

Atom names do have somewhat of a standard format, C1, C2, C3. I am being lazy in the above. If your ligand is present in the RCSB PDB, it will have a chemical component info on the PDB: example and it is often imperative to obey the atom names and atom index from there.
PyMOL or pymol2 module can also add PDB information, therein the residue index is not an integer but a string as it combines w/ insertion code.
In some tools, you lose the bond order, but if you have a SDF and not a PDB, then you are safe. But if you have a reference it's easy to fix.

from rdkit import Chem
from rdkit.Chem import AllChem
template: Chem.Mol = Chem.MolFromSmiles('')
AllChem.AssignBondOrdersFromTemplate(template, mol)

it seem c++ language which i'm not good at, can you provied a python language? the code you provided make error in my pycharm.

It's the compchem module in Python. It runs on PyBind11 C++, so pycharm incorrectly does not like the methods not receiving the instance. #noqa will keep it quiet. I did write the code in the SO editor so I cannot guarantee it runs bug-free.
from typing import List is the only thing missing as far as I can tell.


text:
Generating gene signature for sample classification and survival analysis patient cohort
This is from this paper
The LSC17 score is calculated for each patient as a linear combination of GE of these 17 genes weighted by regression coefficients that were estimated from the training data as follows: LSC17 score = (DNMT3B × 0.0874) + (ZBTB46 × −0.0347) + (NYNRIN × 0.00865) + (ARHGAP22 × −0.0138) + (LAPTM4B × 0.00582) + (MMRN1 × 0.0258) + (DPYSL3 × 0.0284) + (KIAA0125 × 0.0196) + (CDK6 × −0.0704) + (CPXM1 × −0.0258) + (SOCS2 × 0.0271) + (SMIM24 × −0.0226) + (EMP1 × 0.0146) + (NGFRAP1 × 0.0465) + (CD34 × 0.0338) + (AKR1C3 × −0.0402) + (GPR56 × 0.0501). As above- and below-median scores in the training cohort were associated with adverse and favourable cytogenetic risk, respectively, a median threshold was used to discretize scores into high and low groups.
I'm unable to understand or apply the median threshold to create patient group specifically the R implementation here.
In order to address the problem my small dummy data-frame
This is my dummy data
set.seed(123)
nr1 = 20; nr2 = 8; nr3 = 6; nr = nr1 + nr2 + nr3
nc1 = 16; nc2 = 8; nc3 = 10; nc = nc1 + nc2 + nc3
mat = cbind(rbind(matrix(rnorm(nr1*nc1, mean = 1,   sd = 0.5), nr = nr1),
                  matrix(rnorm(nr2*nc1, mean = 0,   sd = 0.5), nr = nr2),
                  matrix(rnorm(nr3*nc1, mean = 0,   sd = 0.5), nr = nr3)),
            rbind(matrix(rnorm(nr1*nc2, mean = 0,   sd = 0.5), nr = nr1),
                  matrix(rnorm(nr2*nc2, mean = 1,   sd = 0.5), nr = nr2),
                  matrix(rnorm(nr3*nc2, mean = 0,   sd = 0.5), nr = nr3)),
            rbind(matrix(rnorm(nr1*nc3, mean = 0.5, sd = 0.5), nr = nr1),
                  matrix(rnorm(nr2*nc3, mean = 0.5, sd = 0.5), nr = nr2),
                  matrix(rnorm(nr3*nc3, mean = 1,   sd = 0.5), nr = nr3))
)
mat = mat[sample(nr, nr), sample(nc, nc)] # random shuffle rows and columns
rownames(mat) = paste0("gene", seq_len(nr))
colnames(mat) = paste0("LSC", seq_len(nc))

They have generated the LSC17 score by doing lasso,here in this case will create a dummy 5 gene score such as this
LSC5 score = (gene1 × 0.0874) + (gene6 × −0.0347) + (gene10 × 0.00865) + (gene22 × −0.0138) + (gene28× 0.00582)

Now to generate or categorize the samples into above or below the median score.
How to implement the above step in R which Im unable to do it?
One suggestion I received is this
 combined into a single variable when you make LSC17. After that, it's just quantile(LSC17,probs=c(0.3333,0.6666))

Which I unsuccessfully followed as I don't know how to implement it in the whole data frame.
Any help or suggestion would be really appreciated
I got an elegant solution which I asked here If there are other ways in R to do this I would love to know


text:
RNAseq: Why would using a more complete scaffold reduce DEGs?
I have a set of RNAseq data for a couple of relatively under-studied species.
They're fungi, of the order Hypocreales. The assemblies we're working with are de novo - reference seqs don't actually exist yet for most of our species. What the collaborator provided was updated genome assemblies in 60-90 fragments as opposed to the original 250-300, but obviously that's still not full resolution. As for the gene families, not a whole lot is known... in fact our species' genes are still not well described and we're mostly working off of predicted orthology.
During the analysis, a collaborator created new DNA scaffolds for these species which we'd now like to use.
However when I redid our RNAseq using the existing reads with the new scaffolds, one of our species had the number of differentially-expressed genes drop by anywhere from 20-50% per condition.
As the reads themselves haven't changed, we weren't expecting significant changes in differential gene expression, and indeed for our other species changes have been minimal.
Does anybody know what might cause this kind of change?
Thank you!

You mention RNAseq reads, but a DNA scaffold. How are you mapping the reads? Could you please show the command line that was used for read mapping, and give some more detail about how you calculated differential expression?

Fungi, ain't a kingdom most of us will know much about. The only thing I can guess at is that in bacteria there has resurgence in methods for de novo assembly specifically aimed at reducing the total number of contigs in an assembly. The 'new wave' is summarised in this response here, Improving prokaryotic assembly with other contig/scaffold-level data?. This is potentially very exciting, except the OP reported they didn't work ;-(
A reasonable guess is your colleague has deployed these approaches to reduce the number of contigs. HOWEVER, what is appropriate for bacteria, may not be appropriate for fungi, alternatively it may be an improvement. Both bateria and fungi are small genomes, but beyond that there could notable differences.
If fungi are dominated by multigene families these might be being removed in the new wave de novo assemblies (they are designed for bacteria). The collapse of multigene families absolutely occurs in other eukaryotes with small genomes, but they have to be really big on multigene families (>10% of genes in their genome). If that was true you would probably know about it.
Alternative, this could simply be poor gene identification in the first assembly resulting in an artificial jump in the number of genes.
Solution? Either way you should be able to figure which genes are different between the assemblies and precisely what type of genes are differing, e.g. are these large multi-gene families. Personally, I would pipe the predicted genes from both assemblies through eggNOG. This might be useful, but not really what it's designed for, but it will return the KEGG nos.  More explicit methods based around comparing the RNAseq reads against the assembly via e.g. bedtools or these guys here recommend  htseq. This might be a better way to rephrase the question.
Thus I recommend:

Finding out which de novo methods were used between first and second fungi genome assemblies. Do any of those methods match this post, Improving prokaryotic assembly with other contig/scaffold-level data?
Consider a more explicit genomics study to examine RNAseq reads against the genome, or comparison between predicted no. genes between genomes.
Having said all that ... once you've done your DEG you  identify the KEGG no. and then simply compare the lists directly to find out which genes are different between the assemblies. This is an easy database operation to perform, e.g. via a merge function (separate question).

I'd check point 1 and then go for point 3 - because it's something I easily know how to do. Thats just my opinion, a purist would likely explore this via point 2.
Conclusion What I am certain about is that for a species where very little might be known it is better to perform an investigation rather than make an assumption that might prove later to be wrong.


text:
Improving conversion of abricate tsv file to gff3 file
Since such a neat solution (abricate tsv to gff3) was provided by Steve, here are few other steps that I am looking to add so that the script progress to logical maturity to be usable by many others.
I have two files - (1) fasta file with .fna extension, and (2) a data file with .tsv extension generated from abricate.
I need to write a gff3 file (containing sequences) by combining these two files (.fna and .tsv)
A .gff3 file generally has three or four components as below. Details here.

Declaration of file type and version
Declaration of sequence region(s) starting from Name of the sequence, followed by single space, then 1 separated by single space, ending in numerical value of the sequence length of the fasta sequence (which is the .fna file).
Sequence features
Fasta sequence (optional)

See below a truncated example
##gff-version 3
##sequence-region NZ_BHEI01000001 1 13890
NZ_BHEI01000001 Prodigal    CDS 29  1588    .   +   0   ID=CEMHHHDG_00001;product=hypothetical protein
NZ_BHEI01000001 Prodigal    CDS 1585    2289    .   +   0   ID=CEMHHHDG_00002;product=Putative beta-lactamase HcpC
....
##FASTA
>NZ_BHEI01000001
AAAAAAAGAATGCTTGCCTACTGGAGTAGTGCCTGTTACTTTTAAACCCACTTTTTTGCG
TTATGCGGATGAGTATTTTTTAAGAGTTGAATTTCAAGATGGAAGTGATGAAATCACTCA
TATAGAGGAGTTGGCAAAATACACAGATCAAA..........
Notice that the sequence names has to match in every instance (declaration of sequence regions, sequence features, and the fasta header at the end)
The .tsv to .gff3 conversion was elegantly solved by Steve.
Now we need to add the declaration of sequence region(s) in the .gff3 file. The declaration comprises of
"##sequence-region"[space]the input file name without the file extension .fna[space]sequence start[space]sequence length[\n]
For the last part, I have used the following awk code on the .fna file to generate the length in desired format.
awk '/^>/ {if (seqlen); seqlen=0;next; } { seqlen += length($0)}END{print "1 " seqlen}' example.fna

This output needs to be appended in the declaration of sequence region line.
Now we need to append the fasta sequence in the following format
##FASTA
>Sequence name as in the parts above in the .gff3 file
AAAAAAAGAATGCTTGCCTACTGGAGTAGTGCCTGTTACTTTTAAACCCACTTTTTTGCG
TTATGCGGATGAGTATTTTTTAAGAGTTGAATTTCAAGATGGAAGTGATGAAATCACTCA
TATAGAGGAGTTGGCAAAATACACAGATCAAA..........
Is there any help in sight, please?

Please use the formatting tools as I requested in your previous question. They really help clarify and that way we can know that the data are correct. It is also much easier than laboriously adding <br> tags.

Thank you @terdon. I do understand that my formatting was not correct. I was also feeling so. Trying to improve! Please bear with for the time being.

Oh hey, that wasn't a rebuke! You spent a lot of effort trying to get your post to look good, it's just that the effort was misplaced because we don't make it very easy to find the right help documentation. I was trying to avoid having you put in so much effort, and ensure that the question is clear and easy to read. So, thank you for the clear effort you put in, even if it did happen to be more than was actually necessary :)

ABRicate can handle multiple input files1 and will set the FILE column in it's output accordingly. Assuming all of the FASTA files in this column exist on disk, you can try the following python solution. It uses pyfaidx, which allows for fast random access to the sequences in the (indexed) FASTA file(s). If the FASTA indexes don't exist, they will be created. It loads only a minimal amount of the input FASTA file(s) into memory. With the following in a file called tsv2gff3.py:
import argparse
import csv
import contextlib
import re
import sys
import textwrap

from dataclasses import dataclass, astuple
from pathlib import Path
from typing import List, Literal, Tuple

from pyfaidx import Fasta

def get_argument_parser():
    parser = argparse.ArgumentParser(add_help=False)
    group = parser.add_argument_group('positional arguments')
    group.add_argument(
        "tsv",
        type=str,
        metavar="FILE",
        help="The ABRicate TSV file",
    )
    group = parser.add_argument_group('output options')
    group.add_argument(
        "-o",
        "--output",
        type=str,
        metavar="FILE",
        default='-',
        help="Write the GFF3 output to FILE (default: stdout)",
    )
    group.add_argument(
        "-f",
        "--force",
        action='store_true',
        help="Overwrite the output file if it exists",
    )
    group = parser.add_argument_group('additional options')
    group.add_argument(
        "-h",
        "--help",
        action="help",
        help="Show this help message and exit",
    )

    return parser

@dataclass
class Feature:

    seq_id: str
    source: str
    type: str
    start: int
    end: int
    score: float
    strand: Literal['-', '+']
    phase: Literal['0', '1', '2']
    attrs: List[Tuple[str, str]]

    @property
    def attributes(self):
        return ';'.join([f'{key}={val}' for key, val in self.attrs])

    def __iter__(self):
        return iter(astuple(self))

    def __str__(self):
        vals = list(map(str, astuple(self)[:-1])) + [self.attributes]
        return '\t'.join(vals) + '\n'

def main():

    parser = get_argument_parser()
    args = parser.parse_args()

    with contextlib.ExitStack() as stack:
        if args.output != '-':
            out = stack.enter_context(open(args.output, 'wt' if args.force else 'xt'))
        else:
            out = sys.stdout

        tsv = stack.enter_context(open(args.tsv, 'r'))
        first_line = tsv.readline().rstrip()
        fieldnames = first_line.lstrip('#').split('\t')
        reader = csv.DictReader(tsv, dialect='excel-tab', fieldnames=fieldnames)

        fasta_files = dict()
        seq_ids = dict()

        out.write('##gff-version 3\n')
        for idx, row in enumerate(reader, 1):

            filename = Path(row['FILE'])
            if filename not in fasta_files:
                fasta_files[str(filename)] = stack.enter_context(Fasta(filename))

            if filename.suffix in ('.fasta', '.fa', '.fna'):
                seq_id = filename.stem
            else:
                seq_id = filename.name

            fasta = fasta_files[str(filename)]
            seq = fasta[row['SEQUENCE']]

            if not seq_id in seq_ids:
                out.write(f'##sequence-region {seq_id} 1 {len(seq)}\n')
                seq_ids[seq_id] = [row['FILE'], row['SEQUENCE']]

            database = row['DATABASE']
            product = re.sub(r'[\[\]\(\)]', '', row['PRODUCT'])
            start, end = int(row['START']), int(row['END'])
            attrs = [
                ('ID', f'{seq_id}_{idx}'),
                ('Product', product),
            ],

            feature = Feature(seq_id, database, "CDS", start, end, '.', '+', '0', *attrs)

            out.write(str(feature))

        out.write('##FASTA\n')
        for seq_id, (filename, seq_name) in seq_ids.items():
            fasta_rec = fasta_files[filename][seq_name]
            seq = textwrap.fill(str(fasta_rec), width=60)

            out.write(f'>{seq_id}\n{seq}\n')

if __name__ == '__main__':
    main()

And pyfaidx installed, the usage is:
$ python tsv2gff3.py -h
usage: tsv2gff3.py [-o FILE] [-f] [-h] FILE

positional arguments:
  FILE                  The ABRicate TSV file

output options:
  -o FILE, --output FILE
                        Write the GFF3 output to FILE (default: stdout)
  -f, --force           Overwrite the output file if it exists

additional options:
  -h, --help            Show this help message and exit

Solution provided by @Steve is elaborate, though import of Literal from typing was causing problem.

@com You might not need pyfaidx if memory is not a concern. It just wasn't clear to me how much memory would be required to load all of the sequences into memory. You can remove the type hints if they're causing problems.


text:
Any existing workflow on Terra that can run multinomial logistic regression?
I need to run multinomial logistic regression on a smaller set of variants e.g. using Trinculo. May I know if there's any existing workflow on Terra that can do this? 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.

What I think is occurring is Terra is the Broad Institute's cloud service. Within that there is a Jupiter notebook and you then just use scikit-learn's (Python) LogisticRegression() using multinomial option. Then is simply ....
from sklearn.linear_model import LogisticRegression
import pandas as pd
model = LogisticRegression(multi_class='multinomial', solver='lbfgs')
df = pd.DataFrame ({'colA':[1,2,3], 'colB':[1,2,5] ... etc ...)
X = df[['colA', 'colB', 'colC']]
y = df[['colD']]
model.fit(X, y)
preds = model.predict(X)
print (preds)

This sort of operation. Does that make sense?


text:
Download all pubmed central article ids that have a keyword appearing in their titles/abstracts and also falling between two dates
I am using the following code to download all pubmed central article ids that have a keyword appearing in their titles/abstracts and also falling between two dates. The code seems to work when I select database = "pubmed" but when I switch to database = "pmc" I get empty list. For example, using the start and end dates as mentioned in the code I get 4 PubMed ids: ['36902449', '36880658', '36583446', '35239250'] and when I check manually, the first one among them (which is '36902449') is also present in PMC database, of course, the id is different (10003631) but same article. The code is expected to fetch that id. I'd appreciate it if someone could point out what I'm missing here.
import os
import datetime
from Bio import Entrez

# Set the email address to identify yourself to NCBI
Entrez.email = "example@gmail.com"

# Set the date range for the search
start_date = datetime.datetime(2023, 3, 1)
end_date = datetime.datetime(2023, 3, 22)

# Format the dates as strings in the correct format
start_date_str = start_date.strftime("%Y/%m/%d")
end_date_str = end_date.strftime("%Y/%m/%d")

# Build the search query with the date range
keyword = "dosage compensation"
# database = "pubmed"
database = "pmc"

# Use esearch to fetch the list of matching article IDs
if database == "pubmed":
    search_term = f'("{start_date_str}"[Date - Publication] : "{end_date_str}"[Date - Publication]) AND {keyword}[Title/Abstract]'
    handle = Entrez.esearch(db=database,
                            term=search_term,
                            retmax=100,
                            api_key = os.environ['ENTREZ_API_KEY'])
elif database == "pmc":
    search_term = f'("{start_date_str}"[Publication Date] : "{end_date_str}"[Publication Date]) AND {keyword}[Title/Abstract] AND open_access[filter]'
    handle = Entrez.esearch(db=database,
                        term=search_term,
                        retmax=1000,
                        api_key = os.environ['ENTREZ_API_KEY'])

record = Entrez.read(handle)
id_list = record['IdList']
handle.close()

# Print the number of matching articles
print(f"Found {len(id_list)} articles published between {start_date_str} and {end_date_str} in {database}")
# print(id_list)

There seems to be a problem with using the [Title/Abstract] field. What you have looks correct and should produce the expected result, since it should capture articles where the search term appears in either the title or in the abstract, or in both. But it does not. Strangely, you get the expected result using {keyword}[Abstract].
A workaround could be to OR the fields with:
f'("{keyword}"[Title] OR "{keyword}"[Abstract])'

Perhaps something like:
search_term = ' AND '.join([
    f'("{start_date_str}"[Publication Date] : "{end_date_str}"[Publication Date])',
    f'("{keyword}"[Title] OR "{keyword}"[Abstract])',
    f'(open_access[filter])',
])


text:
Nextflow Error: failed to read header from "-"
I am trying to run my nextflow pipeline, and have gotten this error:
samtools sort: failed to read header from "-"
I'm not sure why this error is occurring.
#!/usr/bin/env nextflow
nextflow.enable.dsl=2
import groovy.io.*

params.t = 8

process joinFastq {
input:
path(ds_directory)

output:
tuple path("joined_${ds_directory.baseName}.fastq"), val("$ds_directory.baseName")

script:
"""
gunzip -c ${ds_directory}/*/fastq_pass/*.fastq.gz >
 joined_${ds_directory.baseName}.fastq
"""

}
 process minimap2 {
cpus 8
memory "20 GB"
input:
tuple path(fastq), val(baseName)

output:
tuple path("*.bam"), val(baseName)

  script:
   """
  minimap2 -ax splice -uf -k14 -t 8 ${params.genome} ${fastq} > mapped_${baseName}.sam
    samtools view -h -b mapped_${baseName}.sam > ${baseName}.bam
   samtools sort -o ${baseName}.bam > ${baseName}_sorted.bam
   """
   }

I have tried applying the -h flag to samtools sort in the script
   samtools sort -h -o ${baseName}.bam > ${baseName}_sorted.bam

and the same error is produced. I'm not sure why the header flag didn't solve this issue.

This is not a Nextflow error. The problem is that without an input file, samtools sort tries to read from stdin. Using a recent samtools, you can however coordinate sort the SAM and write a sorted BAM using:
samtools sort -o "${baseName}.bam" "mapped_${baseName}.sam"

You may have been intending to pipe the output to samtools sort, which would avoid writing large SAM files and is usually preferable. A minimal example might look like:
params.reads = './*.fastq.gz'
params.reference = './gencode.v43.transcripts.fa'

process minimap2 {

    tag { sample }

    cpus 8
    memory 20.GB

    input:
    tuple val(sample), path(fastq)
    path ref_fasta

    output:
    tuple val(sample), path("${sample}.bam")

    script:
    def cpus = task.cpus > 1 ? task.cpus - 1 : task.cpus

    """
    minimap2 \\
        -ax splice \\
        -uf \\
        -k14 \\
        -t ${cpus} \\
        "${ref_fasta}" \\
        "${fastq}" |
    samtools sort \\
        -o "${sample}.bam"
    """
}

workflow {

    reads = Channel.fromFilePairs( params.reads, size: 1 )
    reference = file( params.reference )

    minimap2( reads, reference )
}

Also, be sure to avoid using absolute paths in your workflow scripts. For your workflow to be portable, processes should only access files in the working directory. The example above avoids this, by using a value channel, to stage the reference FASTA file into the working directory.
In case anyone is wondering, the ternary is to ensure that we do not oversubscribe the node(s) that these jobs land on. Otherwise, each job could use (up to) 900% (800% + 100% = 900%) CPU.


text:
Extract specific nucleotide base position for query and subject from text file
I have an input file looks like this:
> CD630_08510_19
Query  121  TGGCTTGGGTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  121  TGGCTTGGTTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180

Query  241  GCAATTGTGCTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  241  GCAATTGTACTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300

Query  361  CCATACTTTATAGGTATGTTTATTGGTATGGTTGGATAG  399
            ||||| ||||| ||||||||||||||||||| |||||||
Sbjct  361  CCATATTTTAT-GGTATGTTTATTGGTATGGCTGGATAG  398

The input file is the result file after running the command line blast.
It's a txt file generated after blast alignment and in txt file header format are like this.
I am trying to extract the specific nucleotide base for position 129, 249, 366, 372 and 392 from query and subject aligned sequence. python and perl solution will be fine as well. thank you.
For example: I need my output to have format below like this:
Seq_name        base_position       query_base  subject_base    all_subject_base
CD630_08510_19  129                 G           T               TAT-C
CD630_08510_19  249                 G           A
CD630_08510_19  366                 C           T
CD630_08510_19  372                 A           -
CD630_08510_19  392                 T           C

CD630_08510_20  129                 G           G               GATCC
CD630_08510_20  249                 G           A
CD630_08510_20  366                 C           T
CD630_08510_20  372                 A           C
CD630_08510_20  392                 T           C

all_subject_base: are the bases from position 129, 249, 366, 372 and 392 of subject_base column. I want to keep this in one column thats why all_subject_base.
I tried to extract this with grep but have not been able to. I am new in bioinformatics so any suggestion would be very helpful.

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.

Have you tried looking for ways to get alternative output formats from blast that would be easier to parse?

I agree with @MaximilianPress given the approach used. The thing about Biopython is this already done - so it's easy

FYI, the BLAST text format output is for reading, not parsing with a computer program, and is subject to change without notice. I don't recommend using it, especially if you want to do similar things in the future.

Will you always have just one alignment in your file or can you have muiltiple? Does our solution need to deal with both cases?

And why aren't your numbers in sequence? Why do we go from 180 as the last position of the first line to 241 as the first position of the second? Why isn't that 181? Is that actually what your data look like?

Usually, when you run your BLAST alignment you specify the output format (e.g. as described in "outfmt" in https://www.ncbi.nlm.nih.gov/books/NBK279684/table/appendices.T.options_common_to_all_blast/). If, for some reason, you can't re-run the BLAST command (e.g. if this output file was given to you, or if you can't access your query or reference sequence/s) you can parse this 'readable' output format but it's not straightforward.
Here is a potential solution using awk:
cat positions.txt
129
249
366
372
392

cat input.txt
> CD630_08510_19
Query  121  TGGCTTGGGTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  121  TGGCTTGGTTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180

Query  241  GCAATTGTGCTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  241  GCAATTGTACTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300

Query  361  CCATACTTTATAGGTATGTTTATTGGTATGGTTGGATAG  399
            ||||| ||||| ||||||||||||||||||| |||||||
Sbjct  361  CCATATTTTAT-GGTATGTTTATTGGTATGGCTGGATAG  398
> CD630_08510_20
Query  121  TGGCTTGGATTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  121  TGGCTTGGGTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180

Query  241  GCAATTGTGCTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  241  GCAATTGTACTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300

Query  361  CCATACTTTATAGGTATGTTTATTGGTATGGTTGGATAG  399
            ||||| ||||| ||||||||||||||||||| |||||||
Sbjct  361  CCATATTTTAT-GGTATGTTTATTGGTATGGCTGGATAG  398

awk 'BEGIN {
    OFS = "\t"
}

NR == FNR {
    a[$0]
    next
}

$1 == ">" {
b[++c] = $2
}

$1 != ">" {
for (i in a) {
    if ($1 == "Query" && i >= $2 && i <= $4) {
        d[++e] = b[c] OFS i OFS substr($3, (i - $2) + 1, 1)
    } else if ($1 == "Sbjct" && i >= $2 && i <= $4) {
        f[++g] = substr($3, (i - $2) + 1, 1)
        h[c] = h[c] substr($3, (i - $2) + 1, 1)
        }
    }
}

END {
    printf "%s\n", "Seq_name\tbase_position\tquery_base\tsubject_base\tall_subject_base"
    for (j = 1; j <= e; j++) {
        printf "%s\t%s\t%s\n", d[j], f[j], (j % length(a) == 1 ? h[int(j / length(a)) + 1] : "")
    }
}' positions.txt input.txt
Seq_name    base_position   query_base  subject_base    all_subject_base
CD630_08510_19  129 G   T   TA-CT
CD630_08510_19  249 G   A
CD630_08510_19  372 A   -
CD630_08510_19  392 T   C
CD630_08510_19  366 C   T
CD630_08510_20  129 A   G   GA-CT
CD630_08510_20  249 G   A
CD630_08510_20  372 A   -
CD630_08510_20  392 T   C
CD630_08510_20  366 C   T

If you have more positions:
cat positions.txt
129
249
366
372
373
374
375
376
392

awk 'BEGIN {
    OFS = "\t"
}

NR == FNR {
    a[$0]
    next
}

$1 == ">" {
b[++c] = $2
}

$1 != ">" {
for (i in a) {
    if ($1 == "Query" && i >= $2 && i <= $4) {
        d[++e] = b[c] OFS i OFS substr($3, (i - $2) + 1, 1)
    } else if ($1 == "Sbjct" && i >= $2 && i <= $4) {
        f[++g] = substr($3, (i - $2) + 1, 1)
        h[c] = h[c] substr($3, (i - $2) + 1, 1)
        }
    }
}

END {
    printf "%s\n", "Seq_name\tbase_position\tquery_base\tsubject_base\tall_subject_base"
    for (j = 1; j <= e; j++) {
        printf "%s\t%s\t%s\n", d[j], f[j], (j % length(a) == 1 ? h[int(j / length(a)) + 1] : "")
    }
}' positions.txt input.txt
Seq_name    base_position   query_base  subject_base    all_subject_base
CD630_08510_19  129 G   T   TA-GGTACT
CD630_08510_19  249 G   A
CD630_08510_19  372 A   -
CD630_08510_19  373 G   G
CD630_08510_19  374 G   G
CD630_08510_19  375 T   T
CD630_08510_19  376 A   A
CD630_08510_19  392 T   C
CD630_08510_19  366 C   T
CD630_08510_20  129 A   G   GA-GGTACT
CD630_08510_20  249 G   A
CD630_08510_20  372 A   -
CD630_08510_20  373 G   G
CD630_08510_20  374 G   G
CD630_08510_20  375 T   T
CD630_08510_20  376 A   A
CD630_08510_20  392 T   C
CD630_08510_20  366 C   T

NB this won't work if you only have a single position in positions.txt, but if that is a potential use-case you can modify the script to suit, e.g. replace (j % length(a) == 1 ? h[int(j / length(a)) + 1] : "") with f[j].

Jared_mamrot, Thank you so much.


text:
How would you define a set of active genes from a bulk RNA-seq experiment in one biological sample?
I am trying to identify genes that are overall considered to be "actively transcribed" in my cell line of interest for some downstream analysis.
What would the approach be for defining this subset of active genes from an RNA-seq experiments? Any genes with > 1 mapped fragment? Is there a minimum fpkm you would use? or are there some statistics that I can use to identify with some degree of certainty genes that are mostly active in my cell line?
FYI, I am not doing a comparative experiment to identify DEGs, just trying to get a baseline of a group of genes that I can safely considered to be "active" in my cell line.

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

what do you mean by active? To check if a gene is expressed you can look at the counts. The number of counts depends on the biology. Some genes are only lowly expressed because they are toxic at high levels but essential in the cell/tissue. So without more information it is hard to say how to set thresholds. I assume this is bulk RNASeq, did you sequence 1 sample, do you have replicates, biological replicates/technical replicates etc.

Well, that is what I am trying to get at with how people are defining  "active" genes from a  bulk RNA-seq experiment. How would you go about setting a threshold to define a set of MOSTLY active genes in a cell line?


text:
SKAT method in R
I have data that looks like this: 3 column SNPs their gene based on Annovar and a p-value for every SNP. What I would like is to aggregate the p values for every gene.
snps <- data.frame(
  snp_id = c("rs1", "rs2", "rs3", "rs4", "rs5", "rs6", "rs7", "rs8"),
  Gene.refGene_ANNOVAR = c("gene1", "gene1", "gene1", "gene1", "gene2", "gene2", "gene2", "gene2"),
  p.value = c(0.7703884, 0.9648540, 0.9648540, 0.9648540, 0.54, 0.03, 0.03, 0.8)
)

above an example of the data .
I read about the SKAT method -> https://www.hsph.harvard.edu/skat/
and figured it might do the work .
I read about the package here:https://rdrr.io/cran/SKAT/man/SKAT.html
tried to implement it on my data , but got lost as how to perform it correctly:
gene_pvals <- aggregate(p.value ~ Gene.refGene_ANNOVAR, data = df, FUN = function(x) SKAT::SKAT_Null(x)$p.value)

I would be happy if you could share your knowledge in this situation.I don't have information about the correlations between the SNPs but I know that they are correlated, do I have enough data to complete the SKAT method?
also for every SNP i have data on how many patients had the hetro/homozygous encoding(0,1,2).

cross posted: https://www.biostars.org/p/9558525

Unfortunately i didn't recive any  helpful reply in BIOSTARS would be happy for help here


text:
Problem with samtools command with a final hyphen
bwa mem ../01.data/GCF_000146045.2_R64_genomic.fna ../01.data/${sra}_1.fastq.gz ../01.data/${sra}_2.fastq.gz | \
    samtools view -S -b -t ../01.data/GCF_000146045.2_R64_genomic.fna.fai - | \
    samtools sort - >$sra.sorted.bam
samtools index $sra.sorted.bam

I am mapping the downloaded yeast sequences in sra format to the reference genome in a bash shell script.  And I am using samtools to change to bam file.  Then I try to sort the resulted bam file.
I could follow the steps but I could not understand the use of hyphen after the .fai (index file) and the use of this hyphen before the pipe.  Moreover, what is the use of the hyphen after the samtools sort?

could take a look at:  https://www.baeldung.com/linux/dash-in-command-line-parameters

It allows you to chain commands together to avoid creating intermediate files and simplify your workflow:

Samtools is designed to work on a stream. It regards an input file `-'
as the standard input (stdin) and an output file `-' as the standard
output (stdout). Several commands can thus be combined with Unix
pipes. Samtools always output warning and error messages to the
standard error output (stderr).

This is a common idiom in the *nix world, and is actually part of the POSIX specification:

Guideline 13:
For utilities that use operands to represent files to be opened for either reading or writing, the '-' operand should be used to mean only standard input (or standard output when it is clear from context that an output file is being specified) or a file named -.

The hyphen "-" character is commonly used as a placeholder in Unix commands to represent the standard input (stdin) or standard output (stdout) streams. From the POSIX specification:

Guideline 13:
For utilities that use operands to represent files to be opened for
either reading or writing, the '-' operand should be used to mean only
standard input (or standard output when it is clear from context that
an output file is being specified) or a file named -.

In your pipeline, the hyphen after the index file name indicates that the output of the bwa mem command should be piped to samtools view, which will then read it from the standard input stream. The hyphen after the samtools sort command, is not strictly necessary since an input SAM/BAM/CRAM is entirely optional. Note that recent versions of samtools let you coordinate sort the SAM and write a sorted BAM, so you can do away with the call to samtools view entirely. For example:
set -euo pipefail

bwa mem \
    ../01.data/GCF_000146045.2_R64_genomic.fna \
    "../01.data/${sra}_1.fastq.gz" \
    "../01.data/${sra}_2.fastq.gz" |
samtools sort \
    -o "${sra}.sorted.bam"
samtools index \
    "${sra}.sorted.bam"

If you wanted to, you could replace the hyphens with the /dev/stdin special file, which represents the standard input streams. Sometimes this can make more complex commands a bit more readable, but note that not all Unix systems support the /dev/stdin and /dev/stdout special files. Most modern Unix systems do however.
The above uses set -eo pipefail to ensure pipeline errors are not masked. This exits immediately if any command in the pipeline exits with a non-zero status.


text:
Histogram plot using R/python with two variables[BC and Control] and their corresponding gene names
How do we plot histogram in R\Python with two variables[BC and Control] and their corresponding gene names examples is below:

Thank you,

That doesn't look like a histogram but a bar plot. To draw it with ggplot2 and geom_col(), check out this reference: https://ggplot2.tidyverse.org/reference/geom_bar.html

Yeah, I was just given an example through bar plot well, thanks @Cloudberry

How is a bar plot an example for a histogram? They're two different kinds of plots.

Real data is better, particularly if genes names are wanted (complicated simulation) .... anyway there you go. Simulated values under a negative binomial which are pretty realistic to the data supplied (that'll be negative binomial, just not ranked).
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

np.random.seed(2)
s1 = list(np.random.negative_binomial(1, 0.1, 500))
s2 = list(np.random.negative_binomial(1, 0.5, 500))
mydict = {'A':s1, 'B':s2}
df = pd.DataFrame(mydict)
fig, ax = plt.subplots()

heights, bins = np.histogram(df['A'])
heights2, bins2 = np.histogram(df['B'], bins=bins)

width = (bins[1] - bins2[0])/3

ax.bar(bins[:-1], heights, width=width, facecolor='coral')
ax.bar(bins2[:-1]+width, heights2, width=width, facecolor='blue')

Here's horizontal bars (minor differences) with slightly different negative binomial values used ...
# same imports as before

np.random.seed(2)
s1 = list(np.random.negative_binomial(1, 0.1, 500))
s2 = list(np.random.negative_binomial(1, 0.3, 500))
mydict = {'A':s1, 'B':s2}
df = pd.DataFrame(mydict)
fig, ax = plt.subplots()

heights, bins = np.histogram(df['A'])
heights2, bins2 = np.histogram(df['B'], bins=bins)

width = (bins[1] - bins2[0])/3

ax.barh(bins[:-1], heights, width, facecolor='coral')
ax.barh(bins2[:-1]+width, heights2, width, facecolor='blue')
ax.invert_yaxis()
plt.show()

To adapt it you just need to provide an appropriate set of values that will then be fitted to a dataframe.


text:
Running cellranger on scRNASeq data with feature barcoding (10x + antibody capture)
I can't seem to find a clear answer to this question, so here it goes:
I have sequenced scRNASeq + scVDJSeq (TCR) data, which has been sequenced using feature barcoding from 10x genomics, via antibody capture.
There is a very useful cellranger multi command, which allows us to process both this data types in 'one go'.
In order to run it one needs to generate a csv file a priori, such that cellranger is able to identify which files are pertinent to which file type - see this link:https://support.10xgenomics.com/single-cell-vdj/software/pipelines/latest/using/multi
In the link above, with the example given under "Example multi config CSVs", it seems as though there are separate fastq files for the feature barcode, versus for the actual gene expression samples. Whereas, in my case, I only have the standard R1 and R2 fastq files, which is creating some confusion as how to generate this file, and run cellranger.
Now my question is, how does one run cellranger multi with feature barcoding should there only be the standard 2 fastq files? My understanding is that the feature barcode is usually in the R2 read? Is this correct? Would I need to Demultiplex? I am completely lost as to how to proceed. If you know, could you please provide a short example of what your csv files would look like, should you not have separate fastq files for the feature barcode files?
Thanks,

How did you create your fastq files? If you've run cellranger mkfastq with the correct sample sheet, you should have got separate fastq files for each library: RNA-seq, VDJ-seq and feature barcode.

The fastq files were provided by the sequencing company.. We have 4 fastq files (R1,R2,I1,I2) for the GEX & VDJ sequencing, all run via 10x genomics. i.e. I don't have any bcl files, just the fastq's + aligned BAM's

Before running cellranger multi you first need to demultiplex your fastq files to obtain separate files for each library. Unfortunately I've never demultiplexed fastq files so can't recommend any software.

Check your FASTQ headers - if they've been demultiplexed, the index barcodes can be found in the header. The barcodes are usually separated by a + in the comment part (second segment when the header is split by " " blank space) of the header. If you see a match to [ATGC]{10}[+][ATGC]{10} in the header, the first [ATGC]{10} is the I1 barcode and the second is the I2 barcode. You can use FFFFFFFFFF placeholder qual scores to generate "fake" I1 and I2 files and then use that with cellranger.

Thanks for your response. I'm not referring to the I1 & I2 files. I'm referring to feature barcoding technology using antibody capture (or potentially CRISPR). In the examples on the 10x website, they have separate fastq files for these, along with the gene expression fastq's

I see, I apologize for the misunderstanding - I do not know how to help you with feature barcoding.


text:
ucsc-bedgraphtobigwig throwing memory allocation error
I installed the latest ucsc-bedgraphtobigwig (v377) with conda in a new enviroment enviroment inmy WSL2.
conda install -c bioconda ucsc-bedgraphtobigwig
However, when I try using the package, I am getting his error.
free(): invalid pointer
Aborted

I tried downloading the package straight from UCSC, and when I run it, I also get some sort of memory error:
double free or corruption (out)
Aborted

Has anyone experienced something similar? I tried downgrading to v357, and I am still getting the same issue.

Please provide at least the command line you tried and a head of the bedGraph and chromsizes file to make this at least partially reproducible.


text:
Help with analyzing SNP data for three patient groups
I am currently working on a project where I have genetic data for patients classified into three different groups. Specifically, I have information on 26 SNPs for each patient. I have recorded each allele using letters. I am trying to decide on the best approach to analyze this data and would appreciate any guidance or suggestions.
Firstly, I need help formatting the SNP information for analysis. I have recorded each allele using a comma-separated format (e.g., "A,G"). However, other encoding methods may be more appropriate for analyzing SNP data. Could someone please suggest a suitable encoding method for my data?
Additionally, I am interested in performing machine learning and cluster analysis on this dataset to investigate the relationship between genotype and patient group. Please recommend any specific algorithms or packages that may be useful for this type of analysis in either Python or R.
Thank you in advance for your help and guidance.

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

Please describe the input, output and type of ML algorithm you see yourself using for this exercise. From your current text, it looks like machine learning is an appealing concept that you'd like to use in your endeavor, but it doesn't look like you have a concrete idea of what you'd like the machine to learn.

The input would be SNPs information, and the output would be the stage of NAFLD.
Since combining multiple SNPs can help predict the outcome, I am still determining how to analyze these combined correlations. If you have any ideas of statistical analysis I should research, it would greatly help.
About the ML algorithm to use, I would appreciate any recommendations. I could use any supervised ML algorithm but would not know how to find the best-suited one.

Overview In the first instance unsupervised learning is needed - namely PCA - and strongly recommended to identify a training target. The point of this analysis is the resulting clusters of different 'patients' (?).
How to do it? The first part of the analysis is to perform one-hot-key coding.
import pandas as pd
df = pd.read_csv({'A':['A', 'G', 'C'], 'B':['T', 'A', 'G']})
one_hot_encode = pd.get_dummies(df, columns = ['A', 'B'])

This converts A, G, C, T to numerical values and then gets passed to unsupervised learning.
Supervised learning Moving into supervised ML is not trivial, even with a suitable 'target', unless a simple algorithm such as naive Bayes is used. First I'd do the PCA analysis. ...
Edit the training target could simply be the three patient groups.

Thank you for suggesting the one-hot encoding. It has allowed me to perform multiple analyses of the data. As you have said, the training target would be these three patient groups, though I would first like to complete some statistical analysis.
What I truly wanted with this statistical analysis was to find a strong enough correlation between the genetic factors (alone or combined) and the stage of the disease (the three groups) to justify the training of the machine learning model. Still, I need help finding how I can search for this combined correlation.

@IsmaelPrieto Your welcome, although upvotes are much preferred. Firstly bug corrected in the one-hot-key code. Supervised ML classification will identify the genetic factors and assign weights to them. The PCA will assess whether a combined analysis is useful in the ML analysis.

I am sorry, I cannot upvote until 15 reputation. I can, though, accept your answer; that may help.

Thanks @IsmaelPrieto the next question you ask that restriction (20 rep) will (likely) be removed. Basically, in modern ML what you want can singly be done via ML (it's called feature selection/importance). A lot depends on whether you are a pure analyst or you're producing your own data set and analysing that. For pure analytics (its published/public data), the bar is much higher. In your case there will be much more 'leeway'.


text:
Variants list to aminoacids
I have the sequencing data of a Monkeypox virus sample. I have the list of variants extracted from Nextclade and I want to find the corresponding aminoacids substitution.
For example I have the variant C1092T (nucleotide annotation) and I want to automatically convert it into OPG001:D162N (aminoacid annotation).

Please could you provide a reproducible example of the problem you face?

https://www.ensembl.org/info/docs/tools/vep/index.html

Neither C1092T nor OPG001:D162N seem to follow HGVS conventions. You could try to convert them to HGVS nomenclature standards then use modules such as hgvs to get amino acid changes from nucleotide changes. You will need reference sequence identifiers though, so it might be easier to go from chromosomal level changes to AA changes, as cDNA changes depend on the transcript (which may or may not be an issue with viral transcripts).


text:
What is the NCBI's definition of an "atypical genome"?
Using the new NCBI Datasets platform, you can browse the collection of genomes associated with one or more taxa. For example, searching Pseudomonas aeruginosa returns 19,878 genomes as of 29 March 2023.
In the search filtering tab, they give the option to "Exclude atypical genomes". Applying this option to the previous search returns 19,737 genomes, implying that 141 genomes were discarded as atypical.
What criteria are they using to label genomes "atypical"? I've read the documentation for the NCBI Prokaryotic Genome Annotation Process, and my hunch is that atypical genomes are those that are annotated with atypical genes, though I'm not sure.

cross posted: https://www.biostars.org/p/9559023/

Biostars user MirianT_NCBI provided an answer on my crosspost, which I'm copying here:

You can find more information about what constitute an atypical genome in this FAQ. Here's the list of reasons for a genome to be marked as atypical:
Datasets flags genomes with a warning icon and message for the following genome problems:

chimeric - sequences from two different organisms are joined together.
contaminated - sequences from another organism, cloning vectors, linkers, adapters or primers are present in the assembly.
genome length too large - total non-gapped sequence length of the assembly is more than 1.5 times that of the average for the genome in the Assembly resource from the same species, more than 15 Mbp, or is otherwise suspiciously long.
genome length too small - total non-gapped sequence length of the assembly is less than half that of the average for the genomes in the Assembly resource from the same species, less than 300 Kbp, or is otherwise suspiciously short.
hybrid - sequences from a hybrid between different species, strains or isolates.
low quality sequence - long stretches of the sequence have a high proportion of ambiguous bases, are low complexity, or have some other indication that the sequence quality is low.
mixed culture - sequences come from two or more organisms that were not cultured separately.
misassembled - alignment to related genome assemblies or other evidence indicates the assembly is likely to have errors.
partial - the assembly has only partial genome representation.
sequence duplications - assembly has one or more large duplications.
unverified source organism - the origin of the assembly is misidentified.


text:
Using ticks on the x axis or a scale bar when visualisizing a phylogram?
When visuliazing a phylogram is it preferable to show ticks on the x axis or a scale bar?
Given an example tree with branch lengths in newick format (A:0.01,(B:0.02,C:0.05):0.02);, the two options are illustrated in the plot below:

Both will work. The convention is the scale-bar and the rationale is below (which is a bit complicated).
By convention:

the right-hand tree where there is a scale-bar and is the traditional diagram for phylograms.
the left-hand tree scale-bar ("x-axis") is used for molecular clock representations, i.e. not for phylograms, but dendrograms.

In terms of clarity the "x-axis" (or extended scale bar) appears clearer. However, firstly, for a seasoned phylogenecist it makes no difference and simplicity is preferred. Secondly, there's coalescent theory.
The left-hand tree ... the thing is in a molecular clock representation the scale is a reverse axes starting at the right-hand tips. Whether thats time or genetic diversity its a good representation (well for time its essential) because coalescent theory means that there is zero genetic diversity for that taxa at time 0 (the idea is you are extrapolating, or reconstructing, the common ancestor). The representation you presented confuses that convention because 0 starts at the assumed node in the outgroup.
Thus the scale-bar is better.


text:
Manta "--exome: command not found"
I am using Manta (https://github.com/Illumina/manta/blob/master/docs/userGuide/README.md) but whenever I try to specify the --exome command using the following code I get the error below:
The code in my script:
# Accept path to reference fasta file as a parameter:

reference=$1

# Accept input and output directories as parameters:

input_dir=$2

output_dir=$3

# Hard coded Manta Installation Path:

## If someone installs Manta they must specify where the installation is placed:

MANTA_INSTALL_PATH="/path/to/manta-1.6.0.release_src/install"

# Main code:
for file in ${input_dir}/*.bam
do
  # Create subdirectory named after bam file
  sub_dir=${output_dir}/$(basename ${file} .bam)
  mkdir -p ${sub_dir}

  # Use reference fasta path specified by user
  ${MANTA_INSTALL_PATH}/bin/configManta.py \
    --bam $file \
    --referenceFasta ${reference} \
--runDir ${sub_dir}
    --exome
done

# Now run each `runWorkflow.py`

find "$output_dir" -name "runWorkflow.py" -exec {} \;

The runDir is directly below the previous command however, it doesn't seem to write well in this post.
I ran the script via the following command:
./script.sh /path/to/reference.fa /path/to/bam_files /path/to/output_dir

I get the error:
/path/to/my/sample/folder/sample1/runWorkflow.py
/path/to/this/job_scripts/number_by_HPC: line 40: --exome: command not found

How do I get it to run with the --exome command and not get the following error (it works without the --exome command)? I am using exome data so I believe I need to use the --exome command

You're missing a \ in the line before the one containing the --exome, so shell treats it as a new command, which is why you see this error. Add the  \ after ${sub_dir} and your script should not give you this error.
With that correction added in, your code will look like:
# Accept path to reference fasta file as a parameter:

reference=$1

# Accept input and output directories as parameters:

input_dir=$2

output_dir=$3

# Hard coded Manta Installation Path:

## If someone installs Manta they must specify where the installation is placed:

MANTA_INSTALL_PATH="/path/to/manta-1.6.0.release_src/install"

# Main code:
for file in ${input_dir}/*.bam
do
  # Create subdirectory named after bam file
  sub_dir=${output_dir}/$(basename ${file} .bam)
  mkdir -p ${sub_dir}

  # Use reference fasta path specified by user
  ${MANTA_INSTALL_PATH}/bin/configManta.py \
    --bam $file \
    --referenceFasta ${reference} \
--runDir ${sub_dir} \
    --exome
done

# Now run each `runWorkflow.py`

find "$output_dir" -name "runWorkflow.py" -exec {} \;

Thank you Ram. I edited the code above to reflect the change you suggested.

@Indira I think this should have been the accepted answer. Do you still get the error --exome: command not found or something else?

@Indira Did you or did you not run the original command with the backslash? The answer was to your question, not a correction to a typo in your post.

I did run the original commands with a backslash and it did work. Thank you.

In that case, you should have accepted the answer, not changed the question.

I accept the answer, and I changed the code in the question to reflect the correct code that should be used.

That is not how forums work. You let the question have the wrong code and the answer shows the correction. That way, other that made the mistake can figure out the change to be made to their code. I've reversed your edit, accept the answer please.

Oh I will try and keep that in mind for my next post.

@M__ can you clean up the comments please?


text:
Passing strings for password, email authentication and download link
Hi I am downloading the data from a website which requires authentication of email_id and password before the download. They then provide the download link which is used to download dataset. This download link is different every time
### Download Cosmic Files

##First you need to login with your login details and the access key ID may vary with each login ID

##STEP-1
echo "pr@gmail.com:abc123@" | base64
#The above lines gives a authorization key cHIxNEBzYW5

##STEP-2
# get the download link using the authorization key from step 1
#curl -H "Authorization: Basic cHIxNEBzYW5" https://cancer.sanger.ac.uk/cosmic/file_download/GRCh38/cosmic/v97/CosmicSample.tsv.gz

#The above line gives the download link where AWS access key changes everytime like in this case it was https://cog.sanger.ac.uk/cosmic/GRCh38/cosmic/v97/CosmicSample.tsv.gz?AWSAccessKeyId=KRV7P7QR9DL41J9EWGA2&Expires=1679591113&Signature=WPjwDQsMoaLtRA2utiGLNrOSlOA%3D

##STEP-3
# Download file with the link from step 2
#curl -o CosmicSample.tsv "https://cog.sanger.ac.uk/cosmic/GRCh38/cosmic/v97/CosmicSample.tsv.gz?AWSAccessKeyId=KRV7P7QR9DL41J9EWGA2&Expires=1679591113&Signature=WPjwDQsMoaLtRA2utiGLNrOSlOA%3D"

I want to automate the script first by passing used_id and password as an string
I want to replace the AWS key in the step 3 with a string that is generated from step 2

Thank you

The cURL command to obtain the download link returns a JSON containing the URL. This can be parsed on the command line using jq and used to construct a command to download the file. For example, with a list of URLs in a file (one per line) called urls.txt, you could try the following:
while IFS= read -r url ; do
    >&2 echo "Downloading: ${url}"

    curl -sSH @cosmic_header "${url}" |
    jq --arg base "$(basename "${url}")" -r '"curl -o \($base) \(.url | @sh)"' |
    sh
done < urls.txt

Note that the above stores the HTTP header in a file called cosmic_header, which can be created using:
echo "Authorization: Basic $(echo 'username:password' | base64)" > cosmic_header

This avoids exposing the authorization key on the command line when you go to download your list of files, which might be preferable if working in a shared environment.


text:
Analysing the Effect of a Drug on the Morphological Changes of a Cell Type
We came across a project in our lab that no one exactly knows how to approach. Since, I know a little bit of Python programming, this project was assigned to me.
There is a data from a randomised controlled clinical trial with 60 participants. Half (30) are control and the other half are actual patients (case). From each group, biopsy samples were taken before (pre) and after (post) treatment with drug A that targets cell type B in the tissue.
In each biopsy sample, the B cells were detected and classified either as “normal” or “malignant”.
Associated metadata for each biopsy are displayed here (I have only included 12 records for this file):
name    patient_id  arm treatment
111     0       control     pre
112     1       control     pre
113     2       control     pre
121     0       control     post
122     1       control     post
123     2       control     post
211     75      case        pre
212     76      case        pre
213     77      case        pre
221     75      case        post
222     76      case        post
223     77      case        post

●   name: spreadsheet file name
●   patient_id: patient identity number
●   arm: trial arm (‘case’ or ‘control)
●   treatment: treatment condition (‘pre’ or ‘post’ treatment)
Other files (I have only included 20 records for each file: 10 normal and 10 malignant) contain cell detection results from a single biopsy, and each row in a spreadsheet represents an individual detected cell:
●   x: x coordinate
●   y: y coordinate
●   label: classified label (‘normal’ or ‘malignant’)
Just showing an example, file 111 looks like this:
x   y   label
730 724 normal
1962 450 normal
1511 817 normal
1244 455 normal
2529 397 normal
1878 262 normal
2248 369 normal
2007 273 normal
1531 878 normal
1729 834 normal
931  1270 malignant
1282 314 malignant
1630 839 malignant
1543 460 malignant
2493 237 malignant
1311 744 malignant
1999 366 malignant
737 1361 malignant
2252 448 malignant
2620 398 malignant

The rest can be found here, but probably they will not be necessary: https://www.mediafire.com/file/ka7r59kf0swnbnd/OtherFiles.rar/file
I am trying to answer 3 questions here (if you can think of other questions, please let me know):

identifying post-treatment morphological changes due to the effect of drug A.
proposing measures to quantify the changes.
Using appropriate statistical analysis tools to provide insight whether the changes are due to chance or not.

By morphological changes, do you mean the spatial arrangement of cells, or bulk properties of the sample? Since you have coordinates, you could look into "cellular neighborhoods", basically looking at each cell and a small number of its physically nearest neighbors. This could help show if cells are clumping together by type or if they're evenly spread out.

@Nuclear Hoagie I was initially going for spatial arrangment, but at this point i am just trying to come up with the right questions first. Because, i would like to know how to aaproach this problem. I think i should consider both of them.

Yes it's easy you load everything into a dataframe, so each variable is one column and perform supervised machine learning. That will solve all questions.
Hypothesis 1 The target is 'pre' 'post' if there is only one drug in the trial. If there is more than one drug in the trial the drugs would become the target, e.g. if you wanted to know the efficacy between 'drug A' and 'drug B'. What you want to know is the effective of the drug. This would help hypothesis 2.
Hypothesis 2 Estimated cell size is the training target. This is more complicated because I'm not sure what x and y mean: but if x and y are cell measurements then an estimate of cell surface area should be made. If x and y can become a single variate of cell surface area than that can form categorical data ... super size, big, medium large, medium, medium-regular, regular etc ...
This could then become the training target and probably what you are looking for.
The ML feature selection will assess cell size (if 'pre' and 'post' is the training target; or), or 'pre'-'post' (if 'cell size estimation' is the training target). In either case the expected result is it should provide the strongest weight in the analysis.
Finally Finally and importantly, the more variables you can input the better, e.g. the gender of the patient, age of the patient as categorical data (bins), smoking, drinking ... anything you can think of, the more the better. ML should ignore anything that has not impact on 'pre' / 'post' or 'estimated cell size'.
The ML algorithm is more complicated but start with naive Bayes and see how it goes. Beyond this it does start to get complicated however.

To address the comments,

Load everything into a dataframe

A dataframe is a basic unit of operation in Python (pandas) or R. Python also uses numpy arrays.
Everything, means everything - all the data you have, metadata all files. For the calculation the columns involving patient id and patient number, this might be better excluded. I'm not sure - it's complicated because you'll have duplication of records - I would do the calculations with and without this data. You can delete patient id (it's the same as patient number - from a calculation point of view). The x, y is a special case.
Thus all the various tables you describe get loaded into one dataframe, this includes all the patient data such as gender, weight, age, social factors ... it's just like one giant excel table, with rows and columns. The columns are labelled according to the variable e.g. Drug therapy = 'pre' and 'post' data. The dataframe will feed directly to the ML algorithm - a good reason for using it.
x and y
Personally I don't think splitting this data into two columns is a good idea. It will confuse the algorithms - its needs to be one column of data - a single variate. To achieve this will be doable. I would use π to calculate the cell surface area. How that is done for an elongated circle - I don't know - but someone will have done it somewhere, you just need to look it up. Thus 2πr2 if its perfectly round where r = diameter/2 ... elongation must have a derivation its just finding it.
This will absolutely make a difference in your calculations - what I would do is make it to a categorical variable by plotting it and constructing a histogram - otherwise you are into transformations and thats complicated. Categorical variables will also permit logistic regression ML and that can be very useful, because it places each variable into positive or negative weights.

Thanks for your reply. I am not sure if i understand the first line of your response. What do you mean by "load everything into a dataframe, so each variable is one column"? Are you talking about the metadata file or the other files? and which varibales are you talking about?
There are a few things that i should clarify too:
1 - There is only on drug in the experiment.
2 - x and y are coordinates of the each cell detected in the image of tissues in a 2D image.

@EHR answered above below the -----