text:
Can you infer orientation from vcf file data?
I have a dante labs dataset that I would like to cross reference with SNPedia but I am unclear about orientation. I've read the info on SNPedia and other places about it but I just don't get it yet. I realise this is my own stupidity, the explanations feel sound but they seem to be missing some key info that I can't find anywhere.
Anyway without being too tl;dr about what confuses me, I wonder if my dante vcf file has the orientation problem or not because it shows what the reference genome base is next to the sample base. So for example:
1 36571920 rs96067 G A

Chr:1, Pos: 36571920, rsid:rs96067, reference base:G, sample base: A, (another field later in the line defines that both bases are the alt base so this is an (A;A) where the reference had (G;G)).
On SNPedia the genotypes mentioned for rs96067 are C;C, C;T and T;T. Since the reference base in my dante file is a G, which does not show up as a possible reference base on the SNPedia page, does this mean that the dante file is flipped with respect to this gene compared to SNPedia and so I should flip one or other if I want to compare them (flip the base in my file or the data from SNPedia that is)? So I guess flip G->C and A->T meaning my AA becomes TT?
BTW what confuses me about the SNPedia explanation about orientation is this:

Unfortunately if this was instead a SNP where the two alleles are A or T the same flipping logic falls down. We don't (yet) have a way to know for sure if you should be flipped or not, since both forms of your flip rs1234(A;A) and rs1234(T;T) are possible.

Why is AA not the same as TT in this case? Surely this just reflects the direction that the bases were read from, so detection of an AA implicitly means that if we had read from the other direction we would have read TT? I appreciate this is wrong but I don't understand why.
@terdon: Re where I was quoting from, https://www.snpedia.com/index.php/Ambiguous_flip (https://www.snpedia.com/index.php/Orientation also has some info). I've been taking 'flip' as their word for inversion. Re the 'orientation problem' I meant the ambiguity re orientation mentioned in those pages. However you mentioned in your comment that the convention is to always use the + strand. I read up about that and understand now that VCFs use that convention, which answers my question (thanks!). Dbsnp states that it always uses +ve. So like can be compared with like, as long as you use SNP definitions for the correct reference build that your individual SNP VCF has been generated from, and you are using dbsnp as your database. The so called 'flipping' can happen between reference versions/builds and could trip you if you don't use the correct one. Thats what I have gathered to this point, do let me know if its not right.
On SNPedia however many rsid's seem to be defined as minus orientation. They have 2 fields, orientation and stabilised orientation. To quote from the above linked pages:

Orientation indicates the orientation reported in the most current build (currently GRCh38), which is reported below the genotypes in the Reference field. The chromosome number and nucleotide position in that reference build is shown next. StabilizedOrientation is the orientation that is relevant to the genotypes that have been defined in SNPedia for each SNP...

So they are saying that GRCh38 has both plus and minus orientation within and each genotype includes an orientation field? How does this fit with the standard '+ve' convention you mentioned? Maybe these are 2 different things and I just misunderstand. But it would seem far more sensible to simply adjust GRCh38 to all +ve. Is this 'orientation problem' just SNPedias problem? What exactly do they mean on the above linked pages on orientation and ambigous flips?
SNPedia seems pretty useful if it can be made to work - there is lots of associated data on each genotype page, the API is easy to use including for agregated queries etc, but it seems way less reliable than dbsnp despite apparently simply taking their data from dbsnp. Taking the above example I mentioned, rs96067 G>A, dbsnp defines this well, for GRCh37 it shows NC_000001.10:g.36571920G>A. There is also a 'COL8A2 RefSeqGene    NG_016245.2:g.23766C>T' entry which I don't fully understand. However on SNPedia as I said all that shows up are C;C, C;T and T;T genotypes. Maybe they are referring to the COL8A2 C>A entry on dbsnp - or did they scrape this SNP def badly? Or are they showing transcript change rather than genomic as Ram RS suggested?

Try looking up the variant at https://varsome.com. We (I work for the company that builds VarSome, but it is a free tool) integrate data from very many different sources and also show a genome browser which can help put the variant in context. Having a third, orthogonal, source might help clarify the confusion.

SNPedia seems to show transcript base change and not genomic base change, so their C>T is the same as the G>A you note.
I don't understand SNPedia's allele notations completely, other sources are more straightforward. Does their A;A mean A on both strands? Then why do they talk about dbSNP showing C;T when dbSNP would in fact show C>T - they seem to use ; to mean alleles and base change interchangeably, and that is confusing. I recommend sticking to dbSNP for all information not found in SNPedia - something tells me SNPedia is not as trustworthy.

I see what you mean regarding SNPedia, the interchangeable use of the semicolon was actually one of the aspects that was confusing me so I'm glad its not just me. I spent a bit of time checking dbsnp and as you say its always clear. Two problems though, they don't seem to have an API that is as fast easy to use (although I am still trying to get my head round it). And SNPedia summarises related studies which is useful. Re the transcript vs base change, are you sure? I assumed this was just caused by the differences between hg19 and 38 because (once again) its inconsistent (no docs?).

"SNPedia summarizes related studies" - you're going to need to deal with them for that dataset, no way around it. As for dbSNP, you should be able to download a VCF from them that contains all the info you need.

I checked dbSNP and SNPedia for the transcript vs genomic change. SNPedia corresponds to dbSNP's hg38 transcript base change and your input matches dbSNP's hg19 genomic base change. That's the link from your input to SNPedia. I don't know how to proceed with SNPedia though, as their conventions are quite confusing.

Re needing to deal with SNPedia as you say there seems no other way. Although I haven't checked properly it looks like it is a true 'wikipedia' type db with users creating all the summaries and such - theres good and bad in that but when you have so many SNPs to scan I think their approach has value. No need to take what they write as fact but more as an indicator to look more deeply at a particular SNP, perhaps ignore their genotype specific data, instead look at what a particular rsid might have associated with it and if its interesting dig into specific genotypes via another tool.

SNPedia uses ; to separate the two alleles in a genotype. A;A means a (A,A) homozygous genotype. Similarly, C;T is a (C,T) heterozygous genotype. minus indicates the genotype on the GRCh38 strand is (G,A) instead. SNPedia doesn't care about the reference allele.

@user172818 thanks - I don't understand the advantage of this though. Why use C;T and a 'minus' for a genotype that is G;A - why not just use G;A? Can a user of SNPedia always confidently just use the orientation field being minus to flip A-T and G-C in order to get an 'always positive' dataset that will match VCF files and dbsnp without further processing?

From @user172818's comment, it looks like SNPedia lists genotype on cDNA (which makes no sense since transcripts are not double stranded). This semicolon notation is not consistent with 23andMe's online report either. I wonder why this website is popular at all given how different its representation of variants is.

To answer (a small) part of my own question 'Why is AA not the same as TT in this case?' - I believe the answer to this is that the flanking sequence is important here. We can't just flip an AA to a TT in the middle of a strand without having any effect unless we also flip the flanking sequence also. What matters is the relation between the particular base in question and its flanking sequences. I didn't appreciate that before for some reason but its pretty obvious.

SNPedia introduces its own orientation to cut the tie to a reference genome, such that the bulk data is independent of the reference genome in use. The intention is good. The implementation is ok. It is just a bit confusing.
First of all, orientation is always relative; there is no absolute orientation. When you talk about plus or minus, you have to set a "reference" which could be a reference genome, transcription orientation or even something arbitrary.
Then VCF. Orientation in VCF is always relative to the reference that is used to generate the VCF.
In case of SNPedia, "Orientation" is secondary information that is relative to GRCh38, the current human reference genome (PS: this field changes with the reference genome but the internal primary strand system doesn't). rs96067 has two alleles C and T. Because its "Orientation" is "minus", the alleles on the GRCh38 strand is G and A, consistent with the dbSNP information. I checked a few other SNPs. It seems that SNPedia doesn't tell you which is the allele in the reference genome. This is actually fine because what matters is the genotype.
"Stabilized" orientation is relative to an arbitrary orientation stored in the SNPedia database. It is not directly relevant to your question when both you and SNPedia use the same reference genome.
How to use SNPedia? When you have a VCF in the GRCh38 coordinate, traverse rsIDs in the VCF. For each rsID, extract the genotype in the GRCh38 coordinate. Then get the SNPedia "Orientation" and "Geno" based on the rsID. If the orientation is "minus", complement "Geno" (i.e. "A;C" to "T;G", etc). Compare the genoptype from the VCF to the genotype association from SNPedia to predict the effect of genotype.

PS (clarification in response to new comments): if you design a database that ties to GRCh38 from ground up, you may need to flip many genotypes when the reference genome is changed. This process is complex and error prone. With their current design, if the reference genome is changed, they only need to update "Orientation" without changing the core data. This is a thoughtful design, though the UI could be made clearer.

Thanks, very helpful. But I'm unclear about this - 'SNPedia introduces its own orientation to cut the dependency to a reference genome.' vs 'In case of SNPedia, "Orientation" is relative to GRCh38, the current human reference genome.' Those seem to conflict. Surely if orientation is relative to GRCh38 we are still dependant on GRCh38? Otherwise it seems to me + or - is meaningless. What does orientation add? Why not just correct all SNPs to +ve orientation and miss out the orientation field? In fact where did the flip to minus come from - if dbsnp is all +ve, how come snpedia has -ve entries?

@Pete see update. Yes, the two sentences were a bit conflictive. Want to emphasize that "Orientation" is secondary. Their database has its own primary strand system independent of GRCh38.

'Their database has its own primary strand system' - thanks that makes it all fall into place. IIRC even the patch version of the ref dataset matters here (strand can flip from patch version to patch version), so you need your VCF to have been generated against the exact GRCh38-patchX that SNPedia have used to generate their orientation values to be sure you have no errors. How can you map from what they have used (GRCh38-patchX) to GRCh38-patchY or even GRCh37-patchZ? Is the stabilised orientation field relevant there? I can think of a way using dbsnp but is there no simpler way?

Strand flips are common when using different VCF generation software. You've got one of the easier-to-determine variants there, because it's pretty clear that the orientation is incorrect between SNPedia and your VCF file. Illumina is particularly frustrating in this regard because they have a "top"/"bottom" definition, which can sometimes be inconsistent with the "reference positive" convention.
The reason strand for complementary variants can't be [easily] determined is because the usual way to fix strand in a bulk fashion is to check for which of the variant or it's complement are present in the reference database. If the variant is a complementary variant, then both will be found (one as the reference allele, and the other as the alternate variant allele). There are some statistical ways to work round this (e.g. looking at the most common variant in your dataset, and comparing it to the highest allele frequency in the reference database), but those methods will fail if the ancestral background of your study population is substantially different from the reference population (which was the case for the population genetics work I was doing).
The ultimate reference for RS numbers is dbSNP, because that's where the RS numbers are dished out for new SNPs.
Illumina claims that their top/bottom definition allows them to precisely determine complementary strand orientation, but it's useless unless the full sequence context is provided (i.e. the surrounding sequence up until strand can be determined by reverse complementation). Here's the explanation from Illumina:

For SNPs that are not [A/T] or [G/C], A is always on the top strand and T is always the bottom strand. A and T nucleotides are the “A alleles”; G and C nucleotides are the “B alleles”.
If the SNP is [A/T] or [G/C]: Use sequence walking to determine TOP/BOT strands, then assign A/B alleles.

rs903997 ...ATGAGA*A*AGT[C/G]TGA*G*AGTGCA...

By Illumina's definition, the above sequence is designated "top", because if you walk outwards from the variant location, you eventually get to a situation where A is present on one side, and C or G is present on the other side. The reverse complement of this sequence would be designated "bottom":
rs903997 ...TGCACT*C*TCA[C/G]ACT*T*TCTCAT...

This is explained more exhaustively in Illumina's technical note:
https://www.illumina.com/documents/products/technotes/technote_topbot.pdf
When I was puzzling over this, I gave up trying to work this out, and filtered out all complementary variants in my dataset.

I've been really struggling to understand this. Why will both a variant and its complement be found in the ref db if the variant is a complement? Just to check I understand, we are saying the ref is (say) a G, the variant changes that to a T (so G>T). The complement of the T is an A - right? But the ref db still only has GG at this position. Or is a complementary variant a more complex type of variant than a simple base change? Sorry for the probably dumb questions. I'm a beginner...

I've added additional clarification: "one as the reference allele, and the other as the alternate variant allele"

What is this "top/bottom" approach of Illumina's? Can you link to some documentation of it?


text:
RNAseq alignment: best practices for aligning to multiple isoforms?
I have Illumina RNAseq data and would like to maximize my power to find candidate genes that are differentially expressed genes between experimental conditions.
Many of my (de novo assembled and annotated) gene models have multiple isoforms. Because I will only be using "uniquely" aligning read counts (as specified by DESeq2 documentation) I fear that I might be disadvantaging myself by aligning my reads back to the full set of gene annotations.
Would I be wrong to first filter my gff file for only the longest isoform prior to read alignment?

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 program are you using to do read alignment? Please give us the command line (or steps) you used, as the alignment process is different for different programs, and may influence how the transcript models are adjusted.

The fact that an alignment is unique or not is relative to the genome, not its annotation. In other words, a read will be considered a multimapper (and ignored by some programs) if it can align to several genomic positions, but a read will be considered properly mapped if it aligns to a single genomic position, that may or may not correspond to multiple isoforms.
So if you're only interested in gene expression, there shouldn't be any problem with multimappers.
Note that this answer could change depending on what you are aligning on, and how exactly you performed your de novo assembly.

Thanks. I think that you are correct about how "unique" is generally defined. I now feel more confident going ahead while keeping alternate isoforms in by gff file.

If you want to account for changing abundances of isoforms, use a pseudoaligner like kallisto for counting.  And use DEXSeq, not DESeq2, as DESeq2 is intended to be applied to gene counts, not transcript counts.

Yes, and just to be clear, I AM doing read counts per gene. My concern was that those counts could get diluted when alternate transcripts of the same gene are present. However, I believe my understanding of "unique" was wrong and Alexlok's answer below seems correct (from my subsequent searching). Nevertheless, your answer has prompted me to also add DEXseq to my pipeline. So thank you!

Small note, if you do want to look at alternative splicing, in my hands DEXSeq on exon bins was quite inaccurate, I got much better from junction-based tools like spladder and MAJIQ; observation which is essentially confirmed by this recent benchmark. If I'm not mistaken, swbarnes2 is recommending DEXSeq at the transcript-level, which I have no experience or opinion about.

source: Stack Exchange

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

sort: votes


text:
Dealing With Manta Limitations
I am trying to figure out a way to remove the limiting factors that Manta says it cannot handle. To quote from the Manta page:

The following limitations exist on the input BAM or CRAM files provided to Manta:

Alignments cannot have an unknown read sequence (SEQ="*")
Alignments cannot contain the "=" character in the SEQ field.
Alignments cannot use the sequence match/mismatch ("="/"X") CIGAR notation
RG (read group) tags in the alignment records are ignored -- each file will be
treated as representing one sample.
Alignments with basecall quality values greater than 70 are rejected (these are not
supported on the assumption that this indicates an offset error)

Would using a find statement in this way and then making .bai files be adequate to remove the limiting factors that Manta has?:
find ${input_dir} -name "*.bam" -type f -exec sh -c 'output_file="${0%.bam}_filtered.bam"; temp_file=$(mktemp); samtools view -h -q 70 "$0" | sed "/^\(@\|[^[:blank:]]\+[[:blank:]]\+4\)/!s/\t[=*]/\t/g" | samtools view -Sb - > "$temp_file"; mv "$temp_file" "$output_file"; samtools index "$output_file"; rm "$0"' {} \;
# Make the .bai files:

for file in ${input_dir}/*.bam; do
    samtools index "${file}" "${file%.*}.bai"
done

However, I get the error (repeatedly for all my bam files):

[main_samview] fail to read the header from "-".
samtools index: "/path/to/MyBatch1And2/SampleNumber_filtered_filtered_filtered.bam" is in a format that cannot be usefully indexed

I am not sure why I got this error once and then it did not occur again. Why?
Or would this not be enough remove the limiting factors (as there are more to be considered or the code is not written correctly) of Manta?

Assuming you have paired-end Illumina reads mapped using BWA or BWA-MEM (with default parameters), these limitations are usually not much of a concern:

Alignments cannot have an unknown read sequence (SEQ="*")

Best practice for representing alignments in SAM format, suggests that the SEQ and QUAL of secondary alignments should be set to an asterisk to reduce file size. You can either filter these out after the fact, using something like samtools view -F 256 <bam>, or simply avoid writing out any secondary alignments in the first place.

Alignments cannot contain the "=" character in the SEQ field.

An "=" character can be used to denote that the base is identical to the reference base. I've not come across an aligner or downstream process that produces such output. If you know of such a tool or process, please comment below to let me know.

Alignments cannot use the sequence match/mismatch ("="/"X") CIGAR notation

Please see: Why most aligners do not output the "X" CIGAR operation?

RG (read group) tags in the alignment records are ignored -- each file will be treated as representing one sample.

If you have multiple libraries for a single sample, and they differ in some not so small way (perhaps you have a paired-end library and an old mate-pair library, or maybe the distribution of insert sizes/template lengths of each library varies wildly), you may want to consider running each library through separately.

Alignments with basecall quality values greater than 70 are rejected (these are not supported on the assumption that this indicates an offset error)

Short-read sequencing technologies do not produce such high base quality scores. Only long-read technologies, like PacBio and Oxford Nanopore, produce base quality scores above 70 (up to 93).

Ultimately, I don't think you really need to deal with any of these with some command pipeline if you choose your inputs wisely.

As people who don't know how to code and won't choose the inputs as carefully -- will be using my pipeline I cannot afford not to deal with it.

@Indira For a pipeline, I think your best bet would be to start upstream where you get to stipulate the requirements and choose the aligner. The requirements could be one or more pairs of FASTQ files for example. What you want to do is to make sure to fail early. You really do not want to spend several hours or so filtering out all of the long reads from an ONT BAM, only to then have it fall over when the filtered BAM no longer contains any usable data. If you can't start upstream for whatever reason, documenting the requirements might eliminate the majority of issues that arise.

@Indira BTW, I think the problem with your find command is that the sed component is failing. Using single quotes instead of double quotes around the sed expression should solve the problem. Also consider moving the pipeline command inside a separate script file and prepend set -euo pipefail to the start of it. The former should make it a bit more readable, and the latter will catch the errors.

What do you mean start upstream? I am only provided the bams?

I tried to edit the script as well.

@Indira I mean if you can't start with pairs of FASTQ files, documenting the input requirements will go a long way in reducing the number of issues that arise. If you only have BAM files, I would either test to see if the BAM contains suitable alignments or just try to run Manta and include a useful error message in the logs when it fails. People who aren't familiar with your code and don't choose their inputs carefully should still be able to read the error message(s). A notice about the input not being suitable should suffice, and a link to a wiki page could be helpful here.

Regrettably the are not the best at reading error message(s)/will process everything which is why I need to create this script.

@Indira But won't they still get an error when they try to run Manta using a BAM without any read/alignments in it?

What do you mean?

@Indira Well, don't you want Manta to fail when you feed it an empty BAM? If the input BAM contains ONT or PacBio long reads for example, your code will remove all of these.

Considering that there will be batches of samples used I don't want it to fail but to create an out file saying Sample 1 (out of n samples) could not be processed by the main manta script (just not for this limiting factors script which is needed to remove any limiting factors).

@Indira So I think what you want is for Manta to fail, but not your script? In this case, I would skip the pre-processing and just try the command. For example: { >&2 echo "Running Manta for sample x"; ./runWorkflow.py ...; } || { >&2 echo "Could not process sample x"; do_cleanup(); }

Thank you Steve I will try and use that :)

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


text:
Title: Error merging VCF files: "The REF prefixes differ: G vs A"
I'm trying to merge two VCF files, but I'm running into an error message that says The REF prefixes differ: G vs A (1,1) Failed to merge alleles at 1:88143 in sheep2_50KOAR3.1_student_project.vcf.gz. I'm not sure how to fix this error. Can someone help me understand what is causing this error and how to resolve it?
Here is the command I used to merge the files:
bcftools merge sheep1_50KOAR3.1_student_project.vcf sheep2_50KOAR3.1_student_project.vcf -Ov -o sheepcombined_50KOAR3.1_student_project.vcf

Both files have SNP data from the Illumina 50k Ovine SNP, and I believe were generated using the same reference genome. I'm not sure what could be causing this error. Any help would be greatly appreciated.
Thank you!

We can't help you handle data we cannot see. Please [edit] your question and add the two files. We need i) the full headers from both files and ii) the variant(s) that cause the error. Presumably, you have a variant at position 88143 of chr1 that is causing this error.


text:
Parallelize or qsub a bash script
I have a bash script that I would like to parellelize to run on multiple nodes. My goal is to run my python sample_script.py script on pairwise comparisons of samples to see if their variants are a match or not (output is written by echo for the function so we can compare results). I unfortunately can’t post the contents of the Python script, but it is basically comparing germline variants between samples to determine if a sample is a match or not. I’m trying to figure out how to make the script run faster or run multiple pair wise comparisons at once then concatenate results. The samples I'm comparing are listed in the array. Running the script on multiple samples takes too long, so I'd like to qsub the process or run the script on multiple nodes, then concatenate the results. Can someone advise on the best way to tackle this? I don't have much experience with jobs using qsub. Should I consider an SGE array job and if so how would I go about this?
Example code below using bash within jupyter notebook:
%%bash
array=( A037375001 A037646101 A037710601 )
for i in "${array[@]}"
      do
        for j in "${array[@]}"
        do
        python3 sample_script.py --run-sample-id-A "$i" --run-sample-id-B "$j"
        echo run-sample-id-A "$i" run-sample-id-B "$j"
      done
    done > "/outputdir/out.csv"

Thank you in advance!

What's best is quite subjective. My preference would be use Nextflow to abstract away the underlying job submission system and then use the sge executor to submit the jobs to the cluster. The only issue I see, is that it looks like you're feeding in a list of sample IDs, which your script presumably uses to find some VCF files containing the variants. If the script could be refactored to accept file inputs, we can have Nextflow localize them into its working directory. This would ensure your workflow (however trivial) is portable and could be run on AWS Batch, for example. The following example uses the combine operator to create the cartesian product and the collectFile operator to gather the results. With the following in a file called main.nf:
params.vcf_files = './path/to/*.vcf'

process compare_variants {

    input:
    tuple(
        path(sample_A_vcf, stageAs: 'A/*'),
        path(sample_B_vcf, stageAs: 'B/*'),
    )

    output:
    stdout

    """
    your_script.py \\
        --sample-vcf-A "${sample_A_vcf}" \\
        --sample-vcf-B "${sample_B_vcf}"
    """
}

workflow {

     samples = Channel.fromPath( params.vcf_files )
     product = samples.combine(samples)

     compare_variants( product )
         .collectFile( name: 'results.txt' )
         .view { "Saved results to: $it" }
}

And the following in a file called nextflow.config:
process {

  executor = 'sge'
}

And your Python script in a folder called 'bin' in the root directory of your project repository. Nextflow automatically adds this folder to the $PATH in the execution environment. So make it executable, for example using chmod +x your_script.py. Then, with Nextflow installed give it a test:
$ nextflow run main.nf
N E X T F L O W  ~  version 22.10.7
Launching `main.nf` [chaotic_watson] DSL2 - revision: ac224d863b
executor >  local (484)
[bc/a39b27] process > compare_variants (482) [100%] 484 of 484 ✔
Saved results to: /path/to/work/tmp/3e/95330f32edec3f109b197bc836b2ca/results.txt

Thanks Steve! I'm going to try using nextflow. What is stdout referring to? I am getting an error message in my log that this is an unknown variable.

Also, where would I enter in my array info? e.g. the sample IDs?

Steve, I am getting the error message "No such variable: stdout". Any suggestions?

@che625 The above uses the stdout qualifier, but you could also have your script redirect stdout to a file and then use the path qualifier to output the file. You would obtain the sample IDs by either parsing the VCF filenames, or by parsing a CSV containing the sample IDs. For the latter, you could use the splitCsv operator.

Consider using a workflow language such as Snakemake, Nextflow or Cromwell to handle distributing the jobs for you. It will take a little more time investment upfront to set them up, but it will save you many hours and headaches later on.

This is a good suggestion, but as it stands, this is a comment not an answer. How would the OP convert their script to one of your suggestions? What would it look like? Please see Your answer is in another castle: when is an answer not an answer?

You may use asub (array job submitter) for this. Modify your script to myscript.sh:
array=( A037375001 A037646101 A037710601 A037710701 A037753001 A037844701 A038153001 A038172201 A038184401 A038272502 A038298602 A038533901 A038797201 A038944801 A038973401 A039034701 A039048101 A039089601 A039164101 A039303301 A039323801 A039538301 )
for i in "${array[@]}"; do
  for j in "${array[@]}"; do
    echo "python3 sample_script.py --run-sample-id-A $i --run-sample-id-B $j > out-$i-$j.csv"
  done
done

Note that the main difference is this script prints command lines, not executing them. You then run asub with
bash myscript.sh | asub -q myqueue -j myjobname

asub will submit a job array. If you don't like asub, you may submit with
bash myscript.sh | xargs -i echo "echo '{}' | qsub -q myqueue" | sh

This will submit separate jobs, not a job array. It will be slower and less convenient but you don't need to install anything or learn anything.
With either approach, you need to run the following to collate all the results:
cat out-*-*.csv > our.csv

On a side note, this approach can be adapted for GNU parallel to spawn multiple processes on the same machine:
bash myscript.sh | parallel -j 16

GNU parallel can also send processes to multiple remote machines, it isn't only for parallelizing on a single machine. See https://www.gnu.org/software/parallel/parallel_examples.html#example-using-remote-computers

I should clarify that I'm looking to parallelize across samples, so split up running a certain number of samples on different nodes and then concatenate the results. I think parallel will just run the same script on X nodes. Can anyone provide guidance on this?

@che625 I know what you want to achieve. In my answer, parallel executes each line, in parallel, on the same node. I wouldn't recommend to use parallel across node.

So should I divide up my input by positions if I want to use parallel?

@che625 I am not sure what you mean. Anyway, for your original question, just use the "parallel" command line I showed at the end.

Background Just for everyone else SGE = Sun Grid Engine and is the qsub system mentioned in the title of the question.
qsub:
If by parallelising you mean multithreading a single job across loads of cores then the first thing is how many cores per node are there?
Usually it is around 10 (but could be much more). That exact number is super important to know because if you request more cores than exist on a single node the job will never run. At best qsub will refuse the job, at worst it will forever be stuck in the queue. For multi-processing on qsub its ...
qsub -pe omp 8  myscript.sh

Here its requesting 8 cores and myscript.sh will contain the shell command to run the python script. If there's a node with 10-cores 2 are in use, it will then load the job to give its max. capacity.
Thus, qsub does not do course grain-parallelisation, that is MPI. Thus qsub will only parallelise across the cores in a single node. This is without question a limitation for cluster computing.
qsub is more than just submitting to a queue, you need to monitor the queuing system pre- and post- submission. Firstly, this is for the number of existing jobs across that system. Secondly, this is to see whats happening to your stuff.
qstat is the way to understand whats there:
qstat -f # whats what on the cluster
qstat -q

This lists all the queues.
qstat -u username # looks at what a given user is doing
qstat long-queue

Targets availability on a specific queue.
A complete list of qstat is here. Also qusage -l is useful.

"Parallelising" loops
If you are submitting loads of jobs which are all working independently in parallel it better to use a different strategy on qsub. It is better to submit each job separately to the queue via a qsub loop. Thus there would be two scripts the first is the script in the question, the second is a submission loop: this would comprise a qsub argument within a loop that would submit each job sequentially. Thus the array might be better in the qsub submission loop.
Monitoring the queue takes place as per normal via qstat.
Rationale
The reason for this is the way the queue works qsub prioritises single jobs over multi-threaded parallelisation, thus you get your results much quicker.

I should clarify that I'm looking to parallelize across samples, so split up running a certain number of samples on different nodes and then concatenate the results. I think parallel will just run the same script on X nodes. Can anyone provide guidance on this?


text:
How to differentiate DNA fastq and RNA fastq files?
I have 2 sets of fastq files from another collaborator. The first is contains exome data and the second contains  RNAseq data. But both have the RNA in the name but a different ID. How do I differentiate between a fastq file between exome seq and RNAseq?

welcome to the Bioinformatics StackExchange! This is one of those cases where the best solution is to just ask your collaborators. Don't bother with anything technical from the data itself because it's only going to waste your time.

Absolutely agree with James that you should ask your collaborators, but purely as an academic exercise... assuming you have a good reference genome+annotation (e.g. human/mouse), whole exome seq should align to mostly annotated exonic, RNA-seq to exonic+intronic, and whole genome seq to exonic+intronic+intergenic sequences. Sequence duplication levels may also be higher for RNA-Seq with highly expressed genes often 'deeply' sequenced

Adding on this, in WES one sonicates DNA and then captures it. Due to random sonication the coverage around exons should (on the genome browser) look gaussian/bell-shaped. In RNA-seq the transcription start is a precisely defined position, so the coverage at the TSS should bimodally ramp up, same goes for coverage at splice sites, aka it should not be bell-shaped.


text:
How to extract sequences of a specific gene from genomes, with genomes' .fna, .gff and results from KofamScan?
I am trying to do alignment among orthologs and paralogs of a gene but I encountered some problems. I have came up with a workflow, but the sequences I extracted are not correct. They are of hugely varying length, and the alignment is very poor.
I have attached the workflow below, could anyone suggest if I am missing out on anything or if I am interpreting the files wrongly? Really appreciate your help.
My query is: with a list of target genes and bacteria genomes (.fna), for each gene, extract all sequences of it from all genomes and put together into one file for alignment

The answer is very simple, just align from amino acids. Amino acids are not subject to saturation - which is what you are encountering and is smashing up the nucleotide alignment. Amino acid will result in a high quality alignment. Thus you think the sequences are wrong, I think the sequences are simply divergent.
The other issues of variation in length - that depends on the how you've assembled the genes and exactly what project this is. Sinorhizobium is a pretty specific bacteria, so it's pretty clear what the biology is - its quite possibly legume agriculture. However, the precise molecular study is not so clear. Basically, what I'm saying is 'variations' in length are not necessarily 'errors' if this is metagenomics.
You could extract the amino acids from your fna file. Or ... use gff format to get the amino acids via gffread -y protein.fa that sort of thing.
The following might not make much sense ... but hope it helps ... to obtain the DNA alignment ... you convert the amino acid alignment back to a nucleotide alignment by transferring the alignment from amino acids on to the nucleotide data using Emboss' tranalign, amongst many other programs.
The tranalign command is:
tranalign -asequence [nucleotide filename (fasta)] -bsequence [aligned protein filename (fasta)] -outseq

tranalign is also available online, albeit the program is nothing special and its just as easy to write the code.
You need a bit of code to translate the nucleotides to amino acids [I use Biopython translate()] or just use seqotron of MEGAX translate function (if you've a Mac) for a GUI.
Warning If you use this approach the taxa names (strain names) KOs etc ... need to be identical between amino acid and nucleotide files AND the order needs to be identical too. This can be a headache - but if you encounter errors here simply post back as a separate question and I'll forward the code to get around this 'glitch'.
I hope this all sort of makes sense.

HOWEVER, if you want to assess whether KOfamscan is giving the correct results (which is what you are questioning) you could use
KofamKOALA: KEGG Ortholog assignment based on profile HMM and adaptive score threshold available
This is a webserver available here that will do the same thing and should give the same results. The authors consider it an upgrade, but I've not used it.

There are ALOT of other approaches to this problem BTW and all have  over-reliance on KEGG, whilst it is the most comprehensive database available, is not perfect. KEGG per se omits upwards of half of all genes. Nobody else needs to know that but as long as you are aware of this thats fine.

Hi, thank you so much for your detailed answer. Sorry that I missed out another part of the workflow in my original question, but my plan was to translate the .fna to .faa first, do alignment on .faa, and do codon alignment to get the aligned nucleotides. I suspect that the start and end position on the .gff file is not correct, such that I cannot extract the correct nucleotide sequence. I have to get the nucleotide sequence for selection pressure analysis.


text:
Extracting animo acid and nucleotide sequences from KofamScan output and codon alignment
I want to extract the amino acid sequences from KofamScan output, and my workflow is as attached in the picture:
For the analysis I am doing, I need to get the animo acid sequences, align them, and do codon alignment with the corresponding nucleotide sequences, so that I can get aligned nucleotides in the end.
Therefore, for each copy of the gene, I need both the nucleotide sequence and aa sequence. I am extracting these sequences using the .gff file of the genome where the gene can be found, and for each row on the .gff file, there is a distinct combination of the IDs from the .faa and .fna files, and there is start and end position of the target gene that I can use to extract the sequence on the .fna file. (as shown in the 4th step in the workflow)

e.g. NZ_AJQO01000119.1 is the ID from .fna file for nucleotides, and
FIBHGADE03314 is the ID for aa in the .faa file. They are on the same
line in the .gff file, and in that line, that start and end positions
are 60 and 1025 respectively, using the ID NZ_AJQO01000119.1 (as shown in the last step on the workflow).
Therefore, I extract the nucleotide sequences based on the info on
start and end, and do codon alignment with the aligned aa sequences
retrieved using the ID FIBHGADE03314.

However, the extracted nucleotide sequences do not seem correct, while the aa sequences looks good.
I would like to ask if I am interpreting the ID-ID correspondence wrongly? Or am I interpreting the start and end positions wrongly? Are the numbers inclusive/exclusive?
Thank you so much your help.


text:
How to remove orig ident. already data normalized, scaled assigned cluster/
How to remove orig ident. already data normalized, scaled assigned cluster/ Example, I have 4 sample data. like 1-control, 2- treatment with A, 3- treatment with B, and 4- treatment with C,
I want to remove. I want to remove 3- treatment with B, and want to use all 3 ident for downstream analysis. without doing again rescaling and assigning cluster cell type. is it possible?
if yes please tell Sauret r coding so i can run analysis in R

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 [edit] your question and add an example of your input data (use the formatting tools and format it as code by selecting the text and pressing the {} button) and then show what output you would want from that example so we can understand what you need.

The answer depends on how your Seurat object has been created and what the metadata looks like, so please heed terdon's request and provide a representative example of your data and the code you have used.

Can you please add some more filler words into your question? It's very difficult to understand what you're asking in the abbreviated summary format that you've used here.


text:
Where can I find examples of published bioinformatics pipeline packages?
I am a software engineer and I am preparing a presentation to aspiring bioinformatics PhDs on how to use best-practice software engineering when publishing code (such as include documentation, modular design, include tests, ...). In particular my presentation will be focused on publishing "pipeline" code, that is, code that is mainly focused on transforming data to a suitable shape for analysis.
I am trying to find good example of published bioinformatics pipelines that I can point students to, but as I am not a bioinformatician I am struggling to do so, so would like your help. It doesn't matter if the published pipeline is super-niche or not very popular so long as you think it is engineered well so that I can look it up.
Specifically the published code should have: adequate documentation, testing methodology, modular design and be easy to install.

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

@Community I added the list of features that the published package should have to the end of the question.

One quite important point IMO, how are you publishing, e.g. GitHub, Pypy, Docker, Conda? Publishing a data pipeline is not the same as placing monolithic code on GitHub that entirely belongs to the author and has no dependencies subject restricted copyright.

This question is too broad to be useful for the StackExchange format (which works best if there's a single, specific correct answer). I recommend that you ask this at a place where more general discussion is easier, like reddit.

I think the nextflow core (https://nf-co.re) pipelines should satisfy most of your requirements. These pipelines are written in nextflow (https://nextflow.io) which allows scaling from personal to hpc infrastructure rather efficiently.
The core pipelines are usually quite well documented, and have a module design. I have not really looked into the testing methodology, but I think the modules themselves are tested. Running them is easy once nextflow is installed on your system.

The GATK team maintains a GitHub repository with their pipelines: https://github.com/gatk-workflows. They use the Cromwell Workflow Management System, an open source tool created by the Broad Institute.
You are probably most interested in https://github.com/gatk-workflows/gatk4-data-processing, their workflow for aliging fastq files against a reference genome to produce a bam file which can be used for downstream analyses such as variant calling:

gatk4-data-processing
Purpose :
Workflows for processing high-throughput sequencing data for variant discovery with GATK4 and related tools.
processing-for-variant-discovery-gatk4 :
The processing-for-variant-discovery-gatk4 WDL pipeline implements data
pre-processing according to the GATK Best Practices. The workflow
takes as input an unmapped BAM list file (text file containing paths
to unmapped bam files) to perform preprocessing tasks such as mapping,
marking duplicates, and base recalibration. It produces a single BAM
file and its index suitable for variant discovery analysis using tools
such as Haplotypecaller.

If you are starting with FASTQ files visit the seq-format-conversion repository for workflows to convert FASTQs to unmapped BAMS.
The processing-for-variant-discovery-gatk4 provides quick and general processing for sequence data using the latest releases of
GATK. If users are interested in a more elaborate version of this
workflow with quality control tasks and routinely tested for validity
(useful in production environments) then visit the
gatk4-genome-processing-pipeline repository.
The BAM output from processing-for-variant-discovery-gatk4 can be used to perform a variety of other analysis like somatic short variant
discovery, germline short variant discovery, or germline copy number
variant discovery. Visit the GATK Best Practices documentation to
determine what to do next with the BAM files.


text:
DATES: sh: 1: dates_jackknife: not found
I've been trying to run DATES, a software that estimates admixture time, on my Linux Mint machine but I keep running into problems. This is the code to run the program:
./dates -p pardates.txt

The program then runs smoothly, but right before it ends it comes into an error:
Calling fit!...
dates_jackknife -p pardates.txt -l     0.450 -z Assyrian -r 77  -a >Assyrian:log
**sh: 1: dates_jackknife: not found**
##end of dates:       28.806 seconds cpu      627.270 Mbytes in use

I have tried the software on other Linux flavors but it keeps giving the same error, despite the fact I have exported the path to the bin directory every single time. I exported the path with this line:
export PATH=/ubuntu/home/Desktop/DATES/bin

I've also checked my version of perl script (dates jackknife's script) and it's up to date.
I have contacted the authors about the issue but they told me it's an issue with my environment, which I'm pretty sure is not.

it would help to see exactly what the authors said, to diagnose the details.

I suspect the path should not be 'ubuntu/home/Desktop/DATES/bin'. It should be /home/username/DATES/bin, where username would be your ubuntu ID/name


text:
Heritability estimate in LD Score Regression (LDSC)
I am trying to replicate the heritability estimates based on the insomnia GWAS summary statistics using LDSC. However, I have encountered a problem as my estimates seem to be only about half of the original estimates listed in Table S1. Despite my efforts to locate the error, I have been unsuccessful in finding the source of the problem.

Below are the commands I used for female estimates.
zcat insomnia_female_ukb2b_EUR_sumstats_20190311_with_chrX_mac_100.txt.gz |cut -f 6,11,12,15,16,17 -d $'\t'|sed  "1s/STAT/Z/" > female_useful.txt
../munge_sumstats.py --sumstats  female_useful.txt --N-col NMISS --snp RSID_UKB   --a1 A1_UKB  --a2 A2_UKB --out female_insomnia_GWAS_2022 --merge-alleles LDSCORE_w_hm3.snplist --chunksize 500000
../ldsc.py --h2 female_insomnia_GWAS_2022.sumstats.gz --ref-ld-chr eur_w_ld_chr/ --w-ld-chr eur_w_ld_chr/ --out merged_GWAS_results_h2

The GWAS summary statistics were downloaded from
https://ctg.cncr.nl/documents/p1651/insomnia_female_ukb2b_EUR_sumstats_20190311_with_chrX_mac_100.txt.gz
Can anyone tell me where am I wrong? Thanks.


text:
How can I export all and only the publications cited by a specific article?
I want to do a systematic literature review of genetic variant annotation. As part of this I am interested in all the databases available to annotate genetic variants.
I want to export all the references cited by the article with PubMed ID 35943817 (https://pubmed.ncbi.nlm.nih.gov/35943817/), because many of these references are relevant to my review. I can see that they number 93. I want them in a CSV table with their PubMed ID, Title, etc., similar to CSV tables obtained in a search after clicking "Save citations to file", Format "CSV".

Hi @CristianRiccio, your question is fine, but I am not sure this is 'bioinformatics'. The nearest bioinformatics solution is NCBI's Entrez, often used for sequence retrieval, but will retrieve citations from Pubmed - so if you parse the titles in the paper, e.g. Pdfminer you can use those to obtain this data in Entrez, e.g. within Biopython.

You're posting this as a bioinformatics question, which suggests you want code for this. Given that it's only 93 citations, it sounds like you only need to do this once, and they're all in a standard format in the paper. Is there any reason why copying the references from the paper and manually inserting field separators would not be appropriate?

Something like this will work, it did just now for me.  You might want to extract different info from the records though.
This exports to a table (bibliography.tsv) columns containing:

pmid
authors
title
journal
year

max_attempts = 5
def fetch_pubmed_data(pmid):
    Entrez.email = "you@email.com"  # Replace with your email address
    handle = Entrez.efetch(db="pubmed", id=pmid, rettype="xml", retmode="text")
    record = Entrez.read(handle)
    time.sleep(1)
    return record

def extract_cited_articles(pmid):

    record = fetch_pubmed_data(pmid)
    cited_pmids = []

    record_article = record['PubmedArticle']

    if "PubmedData" in record_article[0] and "ReferenceList" in record_article[0]["PubmedData"]:
        for reference in record_article[0]["PubmedData"]["ReferenceList"][0]["Reference"]:
            if "ArticleIdList" in reference and len(reference["ArticleIdList"]) > 0:
                cited_pmids.append(reference["ArticleIdList"][0])

    return cited_pmids

def extract_bibliographic_info(pmid_list, max_attempts):
    bibliographic_info = {}
    bibliographic_info_ls = []

    for pmid in pmid_list:

        authors = title = journal = year = ""
        record = {}

        for attempt in range(max_attempts):
            try:
                record = fetch_pubmed_data(pmid)
                break

            except:
                time.sleep(5)

        if 'PubmedArticle' in record:
            if len(record['PubmedArticle']) >0:
                article = record['PubmedArticle'][0]["MedlineCitation"]["Article"]

                if "AuthorList" in article:
                    authors = ", ".join([f"{author['LastName']} {author['Initials']}" for author in article["AuthorList"] if 'LastName' in author and 'Initials' in author])
                    title = article["ArticleTitle"]
                    journal = article["Journal"]["Title"]
                    year = article["Journal"]["JournalIssue"]["PubDate"]["Year"]

        print (pmid, authors, title, journal, year)
        bibliographic_info_ls.append([pmid, authors, title, journal, year])

    return bibliographic_info_ls

def main():
    pmid = input("Enter the PubMed article ID: ")
    cited_pmids = extract_cited_articles(pmid)
    bibliographic_info_ls = extract_bibliographic_info(cited_pmids, max_attempts)

    df = pd.DataFrame(bibliographic_info_ls, columns=["pmid", "authors", "title", "journal", "year"])
    df.to_csv("bibliography.tsv", sep="\t", index=False)
    return df

if __name__ == "__main__":
    df = main()

Really amazing code. Thank you very much. I just had to change to change cited_pmids.append(reference["ArticleIdList"][0]) to cited_pmids.append(reference["ArticleIdList"][-1]) to get the PMIDs instead of the PMCIDs! Thank you!


text:
GATK g CNV code Aborts Before Running Main Code
I have been struggling to get the code I wrote for running GATK g CNV on my Linux machine (Steps Here: https://gatk.broadinstitute.org/hc/en-us/articles/360035531152 ) to work. Specifically it seems to abort a run before running the main code (I even removed my DetermineGermlineContigPloidy step to no avail). Below is my code and then the error:

## Paths To Be Set

echo "My paths have been set"

### Set Bam Directory:
#### Make sure to set a seperate directory for the batch you are doing

BAM_DIR=/pathToThe/Batch_I_Am_Running

### Set the path to the reference genome:
reference=/pathToMy/primary_assembly.fa

#model_dir =/path/to/my/GATK-gCNV

### Set the path to the singularity image:
SINGULARITY_CACHEDIR=/path/to/my/GATK
singularity_image=$SINGULARITY_CACHEDIR/gatk_latest.sif

### Set Number Of Bins
#### Normally 100 to 200*

NUM_BINS=200

### Setting the number of threads ?? Like Cores ??

NUM_THREADS=20

#______________________________________________________________________________________________________________________________________________________________

## Main code

echo "I am about to start the main code"

for c in "${CHROMOSOME_NAMES[@]}"; do

  # loop over the directories containing the BAM files
  for dir in "${BAM_DIR}"/*/; do

    # get the name of the directory without the path
    dir_name="$(basename "${dir}")"

    # get the name of the BAM file without the extension
    bam_name="$(basename "${dir}"/*.bam .bam)"

    echo "Processing file: $bam_name"

    echo "I am starting Step 1: Preprocessing"

    # Step 1: Preprocessing
    singularity exec $singularity_image PreprocessIntervals \
  -R ${reference} \
      -I ${dir}/${bam_name}.bam \
      -L targets.interval_list \
      --bin-length 0 \
      -imr OVERLAPPING_ONLY \
      -O targets.preprocessed.interval_list

    echo "I am starting Step 2: AnnotateIntervals with GC content"

    # Step 2: AnnotateIntervals with GC content
    singularity exec $singularity_image AnnotateIntervals \
  -L chr${c}.interval_list \
      -I ${dir}/${bam_name}.bam \
      -R ${reference} \
  -imr OVERLAPPING_ONLY \
  -O chr${c}.annotated.tsv

    echo "I am starting Step 3: CollectReadCounts"

    # Step 3: Collect read counts for the samples
    singularity exec $singularity_image CollectReadCounts \
  -I ${dir}/${bam_name}.bam \
  -L targets.preprocessed.interval_list \
  --interval-merging-rule OVERLAPPING_ONLY \
  --format TSV \
  -O ${dir_name}.${bam_name}.counts.tsv

    echo "I am starting Step 4: FilteringIntervals based on GC-content & read depth"

    # Step 4: FilterIntervals based on GC-content and read depth
    singularity exec $singularity_image FilterIntervals \
  -L chr${c}.interval_list \
      --annotated-intervals chr${c}.annotated.tsv \
  -imr OVERLAPPING_ONLY \
  "${dir}"/*.counts.tsv \
      --max-panel-norm 1.5 \
      --min-panel-norm 0.3 \
      -O chr${c}.filtered.interval_list

    echo "Not doing Step 5 DetermineGermlineContigPloidy for now"

    # Step 5: DetermineGermlineContigPloidy (Mandatory Step)
   # singularity exec $singularity_image DetermineGermlineContigPloidy \
 # --model-parameters ${model_dir}/model_parameters.txt \
     # --contig-ploidy-calls ${model_dir}/contig_ploidy_calls.tsv \
 # --contig-ploidy-priors chr20XY_contig_ploidy_priors.tsv \
 # --output-prefix ${model_dir}/ploidy \
     # --standardized-output \
     # $(find "${output_dir}" -name "*.calls.tsv" -printf '%p ') \

#    for file in ${input_dir}/*.tsv
#    do

    # Get sample ID from filename
#      sample=$(basename ${file} .tsv)
#      # Run DetermineGermlineContigPloidy command
#      singularity exec ${singularity_image} DetermineGermlineContigPloidy \
#        --input ${file} \
#        --contig-ploidy-priors ${contig_ploidy_priors} \
#        --output-prefix ${output_dir}/${sample}_ploidy \
#        --standardized-output
#    done

    echo "Creating shards via bins"

    # Creating shards via bins
    singularity exec $singularity_image IntervalListTools \
  SplitIntervals \
  -R ${reference} \
      -L chr${c}.interval_list \
  -scatter ${NUM_BINS} \
      -o ${dir}/${c}_shard.interval_list

    echo "Starting the Shard Loop that will run GermlineCNVCaller and PostprocessGermlineCNVCalls"

    # Loop over the shards and run GermlineCNVCaller and PostprocessGermlineCNVCalls
    for shard_file in "${dir}/${c}"*_shard.interval_list; do

      echo "Starting Step 6: Run GermlineCNVCaller (Mandatory Step)"

      # Step 6: Run GermlineCNVCaller (Mandatory Step)
      singularity exec $singularity_image GermlineCNVCaller \
    --run-mode COHORT \
    -L "${shard_file}" \
        -I $(find "${dir}" -name "*.tsv" -printf '%p ') \
        --minimum-interval-median-percentile 5.0 \
        --maximum-zeros-in-sample-percentage 40.0 \
        --annotation-cutoffs "QD=2.0" "SOR=2.0" "FS=10.0" "MQ=50.0" \
        --interval-merging-rule OVERLAPPING_ONLY \
        --output-prefix $(basename "${shard_file}" .interval_list)

      echo "Starting Step 7: Run PostprocessGermlineCNVCalls (Mandatory Step)"

      # Step 7: Run PostprocessGermlineCNVCalls (Mandatory Step)
      singularity exec $singularity_image PostprocessGermlineCNVCalls \
    --model-shard-path "${shard_file}" \
        --cnv-call-path $(basename "${shard_file}" .interval_list).calls.tsv \
        --output-dir postprocess-sm \
        --scatter-count ${NUM_THREADS} \
        --num-scatter-plot-columns 8 \
        --max-zeros-in-sample-percentage 40.0 \
        --annotation-cutoffs "QD=2.0" "SOR=2.0" "FS=10.0" "MQ=50.0" \
        --interval-merging-rule OVERLAPPING_ONLY \
        --minimum-contig-length 46709983 \
        --minimum-interval-median-percentile 5.0

    done

  done

done

    #______________________________________________________________________________________________________________________________________________________________

The ".out" file I received after the run aborted is:

My paths have been set
I am about to start the main code

Why would it abort so soon (I changed the number of Gs used and that did not have any effect)?

Please don't dump 168 lines of code on us and expect us to debug it. You need to make this more focused: give us a minimal script that still produces the error. You haven't even put in the minimal effort to remove commented lines of code here! In this case, the error is trivial, but as a general rule, you need to help us help you.

I will try. I am still getting used to coding. Shouldn't commented lines remain?

No worries, we all start somewhere! Just remember the first rule of debugging: print all the things! By far the most common source of bugs is that the developer assumes their variables have one value, but actually the variables have another value. So the first step, always, should be to print out the variables, and to add simple print statements like echo "I got to here!" at different places in your code so you can figure out exactly where the problem occurs.

That is what I was trying to do with the echo statements.

You are using ${CHROMOSOME_NAMES[@]} but you never define that, so when your code gets here:
for c in "${CHROMOSOME_NAMES[@]}"; do
...
done

It doesn't enter the for loop because "${CHROMOSOME_NAMES[@]}" is an empty array. It's the same thing as doing this in a shell session:
$ for c in "${CHROMOSOME_NAMES[@]}"; do echo "This is chromosome $c"; done
$

Just like your script, that will produce no output since there is nothing for it to loop over.

Thank you for this I will try it out seems that I removed my code CHROMOSOME_NAMES=( $(seq 1 22) X Y)


text:
How to calculate the Maximum Ladder Distance (MLD) of ssRNA from Dot-Bracket notation?
The ladder distance (LD) between residue i and residue j is the minimum number of base pairs (created due to secondary structure of the ssRNA) one has to cross when going from i to j.
The Maximum Ladder Distance of a ssRNA molecule is the maximum ladder distance between any two of the residues of this molecule[*].
Given the secondary structure of a ssRNA molecule in dot-bracket notation (e.g. by RNAfold), how can I calculate the Maximum Ladder Distance of this molecule?
Below is an illustration from [*]. The MLD is 11, since this is the number of basepairs on the yellow line.
Edit: The dot-bracket notation for this molecule would be
..(((((.......))...((.((((((.......))))))..)).)).

*: Yoffe AM, Prinsen P, Gopal A, Knobler CM, Gelbart WM, et al. (2008)
Predicting the sizes of large RNA molecules. Proc Natl Acad Sci U S A 105:
16153–16158

Please [edit] your question and post the dot-bracket representation of the above structure. Use the code formatting button ({}) to format it in a monospace font for clarity.

There are different implementations for calcuating the $\rm MLD$. One possible approach is to consider $\rm MLD$ as the heighest value of all possible mountain plots, generated by shifting the beginning of your sequence. The highest value of a mountain plot tells you the maximal ladder distance that can be reached from the first nucleotide. But the longest path might not begin at the first nucleotide, and that is why you have to go through all possible shifts of the sequence.
Here is an example of this implementation in Python. The algorithm is as following:

Convert dot-bracket notation to a list of base pairs (dot2pairs).
Calculate the first mountain plot.
Calculate required shifts. This is not strictly neccessary because you can swipe through all possible shifts = list(range(0, length)), but it gives better performance.
Iterate through all shifts, calculate the MLD at each step, and return the maximal one.

import numpy as np
from collections import deque

def dot2pairs(dot):
    stack, pairs = deque(), []

    for i, n in enumerate(dot):
        if n == '(':
            stack.append(i)
        elif n == ')':
            pairs.append((stack.pop(), i))

    return np.array(pairs)

def dot2MLD(dot):
    # Convert to a list of base pairs
    pairs0 = dot2pairs(dot)

    if len(pairs0) == 0:
        return 0

    MLD = 0
    shifts = [0]
    length = len(dot)
    seq = np.zeros(length, dtype=int)
    curr_mount = None

    while shifts:
        # Shift the sequence
        pairs = (pairs0 - shifts.pop()) % length

        # Find positions of opening and closing nucleotides
        opens, closes = np.min(pairs, axis=1), np.max(pairs, axis=1)

        seq.fill(0)
        seq[opens] = 1
        seq[closes] = -1

        mountain = np.cumsum(seq)

        MLD = max(MLD, int(np.max(mountain)))

        # In first iteration, calculate all required shifts
        if curr_mount is None:
            curr_mount = mountain[0]
            for i, m in enumerate(mountain):
                if m == curr_mount:
                    shifts.append(i)
                else:
                    curr_mount = m

    return MLD

Let's test it on your structure:
dot2MLD("..(((((.......)))...((.((((((.......))))))..)).)).")

11

And on a larger structure:
dot2MLD(s)

1141

It took ~ 5 seconds for this 17481 nt long structure (and it takes ~ 20 seconds if you don't precalculate requires shifts).


text:
How to get enriched pathways in the data using continous statistic measure?
I was doing pathway enrichment analysis using the below code
p_value = scipy.stats.hypergeom.sf(len(intersection) - 1, geneset_size, len(pathway), len(hits))
        if p_value < 0.05:
            for gene in list(pathway):
                results = results.append({"Gene": gene, "Model": model, "Enriched Pathway": Pathway_name}, ignore_index=True)

where,
hits = df["Gene"][np.logical_and(df["Model"] == model, df["p_value_anova"] < 0.05)]
intersection = pathway.intersection(hits)
Model = model name (like A, B, C) in input file
Gene = Gene from the input gene
Enriched pathway=Pathways that are enriched from the Pathway input file provided

Now the problem is that earlier this checking pathway enrichment on the basis of p_value in hits was working but now using cutoff (0.05) is resulting in loss of all pathways. I am looking to test the enrichment using the statistical test that will work on continous measure than binarizing the P-value.

Hi @Megha, do you have a small amount of sample data?

Just to mention I personally think there is a higher loop above the loop stated above.
I've not looked at each bit of code in detail but
there's a bug is on line 4:
results = results.append({"Gene": gene, "Model": model, "Enriched Pathway": Pathway_name}, ignore_index=True)

Firstly, this is fashionable coding because you are dumping straight into a pandas dataframe (described below). The correction is simply ... declare pandas in advance
import pandas as pd
results = pd.DataFrame()
p_value = # code missing
if p_value < 0.05:
    for gene in list(pathway):
        results.append({"Gene": gene, "Model": model, "Enriched Pathway": Pathway_name}, ignore_index=True)

The thing is if there is a higher loop above list(pathway) then results = results.append(.... will over-ride the last set of results. This would explain why results 'disappear' because if the last data in the loop is a blank - or very little output, all previous stuff gets deleted.
I personally would use a dictionary here. Thus
# declare the dictionary before all loops
results = {"Gene": [], "Model": [], "Enriched Pathway": []}
p_value = # code missing
if p_value < 0.05:
    for gene in list(pathway):
       genedict = {"Gene": gene, "Model": model, "Enriched Pathway": Pathway_name}
       for k,v in mydicti.values():
           results[v].append(v)

... you could then transfer this to pandas once the loops are complete via,
df = pd.DataFrame(genedict)

I do suspect there's other issues, but as a general rule results = results.append(... is better declaring results first and then appending to results, unless you've a very specific reason (e.g. writing each loop into a separate file). The pandas thing is kinda clever for big data sets, because pandas is fast, whilst dictionaries are slow - for big data (like >>500 000 big). The interchangeability between dictionaries and pandas means it's kinda easier conceptually (and code-wise) to switch at the coders convenience in my personal opinion.


text:
Test for differences between groups of samples
Sorry if the answer to this should be obvious.
I have RNA-expression results from 24 samples which can be divided into 6 groups, (wildtype and two different mutants at two different ages) with a total of 24 samples. I have made a PCA plot and when describing it in my thesis said separation between samples is minimal. One of my examiners has said

You state that the separation between samples in the principal component analysis is minimal, this needs to be quantified to differentiate between minimal but statistically significant, and not significantly different.

I have no idea how he expects me to test for significance. I know how to test for significant differences between individual samples, and for differences per gene/principle component. I don't know how to group samples together while testing for significance against multiple things or if I should use genes or principle components (probably the latter).
I'm writing it in R so if anyone knows which R function can do the kinds of statistical tests I need that would be a huge help. Alternatively is this something which I shouldn't really do and I should push back on?
I'm already performing significance testing using DESeq2. To the best of my knowledge that only tests for whether a single gene is significant. It doesn't test for whether the samples are significantly different as a whole which is what I'm after.

Hi Sethzard I think you need to do is present the experimental design. The problem at present is it's not clear what the criticisms are because we don't know your experiment for "grouping samples together". @RamRS is correct you cannot use PCA. At least you cannot use traditional PCA as used in RNAseq. It's also clear "I don't know how to" is this code or test?

PCA does not do any testing and has no p-values. You may state that a separation of less than m units in PC > i is not significant while a separation of greater than n units in PC < i is significant as PCs < i explain X% of the variation (i being 2 or 3). That's how you can back your statement on the separation being minimal.
That is, in this sample plot:

PC2 distance of < 0.05 might be called trivial compared to PC1 distance > 0.5 if PC1 explains 80% of the variance. This is just 2 PCs though, so your actual case will need to be more robust.

I know that PCAs don't do any testing on their own. I wasn't sure if they could be used as part of a statistical test which does test for significance or if I'd need to use the raw values. Thank you for the suggestion of other ways to do it.


text:
Combine separate fasta files from a list.txt
I have about 10 fasta files that I need to pull from a giant folder of other fasta files (has about 100 files). I have a list of the files I need in a .txt file separated by line and I want to pull out the 10 fasta files listed and combine them into one big file.
I tried this command at first:
xargs < names.txt cat > names.fasta
But it gave me this error:
cat: 'read_1.fasta'$'\r': No such file or directory
For all the files except the last read in the list which it managed to add to the names.fasta .
I'm not sure what the error message means given I can find the files in my folder so they should be there. I can do this manually but I know it'll take a while to do.
EDIT: I changed the list to be separated by space and it actually pulled out all the files I needed! However, is there a way to account for reading in the names by newline instead of space??

Please always mention what operating system you are using. Your error suggests you are using Windows and have Windows line endings (\r\n) and that is essential information when trying to manipulate files and parse text.

This should work
xargs < names.txt cat > names.fasta

Or else try,
cat names.txt | xargs cat > newdir/names.fasta

What I think is within names.txt there's a typo or a missing file and that's what it's complaining about.
Perhaps try using
ls > tmp
cat names.txt | sort > tmp2
diff tmp tmp2

If there are any differences that should pick it up.

Just to address an edge case, run dos2unix names.txt first and then a modified version of @M__'s command
cat names.txt | xargs -I v_file cat v_file >> names.fasta

Based on your edit, it looks like your working platform has different line endings compared to the platform on which names.txt was last edited/saved.


text:
how to handle variable output from nextflow
I have this output which is going to do bwa next which is from join operator
[DF122_SO11378_A03, DF122_SO11378, DF122_SO11378_A03.norm.merged.fq, DF122_SO11378_A03.final.gfa, DF122_SO11378_A03.final.fasta]

in this case, pair_id is DF122_SO11378_A03,
plate_id is DF122_SO11378, fq is DF122_SO11378_A03.norm.merged.fq,
DF122_SO11378_A03.final.gfa  is gfa
DF122_SO11378_A03.final.fasta  is fa
I have following code like this
process bwa {

     input:
     tuple val(pair_id), val(plate_id), path(fq), path(gfa), path(fa) from align_in

     output:
     path "*bam" into bam
     tuple val(pair_id), val(plate_id), path ("*.csv") into metrics

     """
     if [ -s $fa ]; then
     bwa index $fa

     bwa mem $fa $fq \
      | samtools view -bh -F2048 - \
      | samtools sort > ${pair_id}.bam

     samtools depth -a ${pair_id}.bam > ${pair_id}.depth.csv
     samtools view -c ${pair_id}.bam >${pair_id}.count.csv

     else
     touch ${pair_id}.depth.csv
     touch ${pair_id}.count.csv
     touch ${pair_id}.bam
     fi
     """
}

However, now my upstream work flow changed a little, for some sample they have more than one gfa and fasta output like this
[DF122_SO11378_A04, DF122_SO11378, DF122_SO11378_A04.norm.merged.fq, [DF122_SO11378_A04.2.final.gfa, DF122_SO11378_A04.final.gfa], [DF122_SO11378_A04.2.final.fasta, DF122_SO11378_A04.final.fasta]]

in this case, pair_id is DF122_SO11378_A04,
plate_id is DF122_SO11378, fq is DF122_SO11378_A04.norm.merged.fq,
[DF122_SO11378_A04.2.final.gfa, DF122_SO11378_A04.final.gfa] is gfa,
[DF122_SO11378_A04.2.final.fasta, DF122_SO11378_A04.final.fasta] is fa
But some samples still have one gfa and fasta as shown earlier.
How can I adapt my bwa process to handle such situations?
I am using DSL1.

I think I can use transpose operator to make the list flat, and use the downstream pipeline. So the problem is currently solved.

You can put the details in an answer, and then accept your own answer as the solution.  That'll help anyone in the same situation later.

Note that there's no reference to the GFA file in the script block of the process. Also, testing if the FASTA file is empty is probably unnecessary. Ideally, this should be performed upstream so as to not execute tasks that just create empty (placeholder) files.

If you're looking for the transposition, you could use the transpose operator to get what you want. But to handle the variable output, you would need to first check what types you have. For this, you could use something like the following for example:
input_ch
    .map { pair_id, plate_id, fastq, gfa, fasta ->

        def gfa_files = gfa instanceof Path ? [gfa] : gfa
        def fasta_files = fasta instanceof Path ? [fasta] : fasta

        tuple( pair_id, plate_id, fastq, gfa_files, fasta_files )
    }
    .transpose()
    .view()

Note that the old DSL1 is no longer supported.


text:
When it is stated in your consumer genetics testing report that you have x% (lets say 1.3 percent...) Neanderthal DNA, what is actually meant?
This question has also been asked on Biology SE
When a raw number like 1.3 percent is stated, it is confusing without context,and can be misinterpreted. So does 1.3 percent mean that 1.3 percent of the letters/nucleotides/bases of my genome consists of established (likely/high confidence) Neanderthal markers or SNPS? Does it mean 1.3 percent of my entire genome consists of these Neanderthal markers, i.e 1.3 percent of 3.1 billion or 6.2 billion base pairs, or is it 1.3 percent of protein coding regions, or exons?**
How do they obtain these markers for Neanderthal ancestry to begin with? I know the Neanderthal genome has been sequenced from bone fragments. Do they compare SNPs and variable regions in the genomes between Neanderthal and Modern Human population genomes to see the relative frequencies of these variations/SNPs. Then likely Neanderthal markers are flagged by statistical models and algorithms..for example combinations of SNPs that appear infrequently or hardly at all in the reference Sapien population genome but occur frequently in the Neanderthal reference genome..that marks these SNPs/variable regions as of likely Neanderthal origin? Could you please describe some of these statistical methods and algorithms? Can anything else explain the "markers", or SNPs other than interbreeding with Neanderthal (such as mutations)? Please illustrate using some hypothetical alignments like AAGGGT..vs GAGGGT..Please give a clear, concise and authoritative answer of high confidence. It is okay to be detailed and technical, although building up to complexity so that it can be understood with relatively little familiarity with population genetics and Mathematics.

I ran a basic search of 23andme and it looks like they use known variant loci from the Neanderthal genome. Also, they expect a 2-2.5% overlap. See:

https://customercare.23andme.com/hc/en-us/articles/212873707-Neanderthal-Ancestry-Report-Basics
https://blog.23andme.com/articles/new-neanderthal-report (Section "The Science behind the report")
https://permalinks.23andme.com/pdf/23-05_neanderthal_ancestry_inference.pdf (Eric Durand's white paper on the technical details)

So the percentage they give corresponds to the percentage of the specific set of variant loci that a test sample has in common with a Neanderthal reference genome? What does the actual percentage mean, exactly?

@terdon I believe the percentage gives the fraction of your genome derived from a Neanderthal genome.

@user172818 yeah, but fraction of what, exactly? Percent of similarity between coding regions? Between known transcript regions? Genes? Across the whole genome including non coding and intergenic areas? That is what the OP is asking, as far as I can tell, and the answer isn't clear (to me) so far. Perhaps the answers are in the given links, but if so, they are not here.

@terdon Whole genome

Do you mean the results from the 'whole-genome' SNP (genotyping) chips?

Based on reading a similar post it seems like if the report says I am 1.3 percent Neanderthal  then it  does NOT  mean 1.3% of my  base pairs or SNPs of my genome (or exons or protein coding region or total number of my SNPs) has a "Neanderthal" label inferred, rather some kind of statistical inference is made as a whole.
The  author suggests when it comes to 23AndMe  that my genome is segmented into chunks of a chromosome (containing a set of SNPs, 23AndMe does genotyping only reads parts of chromosome where bases are known to vary.)
Technical details could be referred to the following sites:  https://academic.oup.com/bib/article/19/5/765/3038466
https://www.cd-genomics.com/the-principles-and-workflow-of-snp-microarray.html
https://www.illumina.com/techniques/popular-applications/genotyping/snp-snv-genotyping.html
Since 23andMe is using v5 genotyping chip, the fraction of Neanderthal variant could be calculated.  Moreover, the basis of v5 genotyping chip is related to the concept that DNA-based polymorphisms (SNPs or SNVs) which differ substantially in frequency among human populations can be employed to infer ancestry.  With 23andMe, an average of 236 Neanderthal SNPs/SNVs could be obtained with the analysis based on a selected group of consented customers.
And these chunks are fed into a Machine learning classification  algorithm like Support vector machine (SVM). The Machine learning algorithm was trained chunks of genome of  reference population sets with chunks of known ancestry (like "Scandinavian" or "South European" or "Sub Saharan African").
Then the classification algorithm makes a multilabel prediction... lets say the algorithm predicts that a particular  segment is 95 percent likely to be human and 5 percent Neanderthal..and maybe the next segment has a  2% probability of being human and 98 percent Neanderthal..Then the algorithm looks for segments of high confidence (>90 percent or 70 percent whatever threshold you set) and then add them up.
For example if in a toy example genome has 9 segments of the same size has a  greater than 90 percent probability of being estimated "human", and 1 segment with greater than 90 percent of being estimated  Neanderthal by the classification algorithm than I would be  considered to be 10 percent Neanderthal...
The author suggests this is the method used to infer ancestry like "Norwegian" and "South European" and "Sub Saharan African" by 23And Me..but will it be similar for inferring Neanderthal ancestry?
Calculating ancestry percentage is not seem as  simple as counting up number of "Neanderthal" SNPs/Markers and dividing by your total SNPs or base pairs in your genome/exon...Rather an inference of ancestry mix needs to be made based on reading SNPs on segments of chromosomes and comparing it to reference data sets via classification algorithms.

Also could someone clarify what "segments" or chunks mean. when someone mentions "..take the genome of target individuals and split it up into chunks along the genome". I know 23 and Me does genotyping not WGS or whole exome sequencing..so these segments or chunks must consist of a a set of SNPs..presumably, rather than a continuous segment of the Genome including invariable/conserved region...?

You are overthinking it. What method they use doesn't matter much. The goal of ancestry estimate is to get genome fraction.

'The most commonly used criterion for SNP (genotyping) chip evaluation is global coverage, defined as the fraction of common SNPs that are tagged by the SNPs on the chip.'  Do you mean the fraction of common SNPs is the same as the number of SNPs tagged by SNVs on the SNP (genotyping) array?

@SanjayBiswas  You are right.  Current 23andme is using their v5 genotyping array and genotyping arrays could measure up to millions of SNPs across the genome as illustrated in this table --> https://www.nature.com/articles/s41431-021-00917-7/tables/1

@SanjayBiswas  Please also refer to this slide shown in the powerpoint of Lynn Jorde --> https://youtu.be/PzG99bexziA?t=3712

How do they obtain these markers for Neanderthal ancestry to begin with? I know the Neanderthal genome has been sequenced from bone fragments

Since ancient nuclear DNA is scarce, enough preserved molecules are needed.
The level of scarcity about ancient nuclear DNA : often, 1% or fewer of the total DNA fragments from an extraction among soil microbes and other sources
After extraction and sequencing, distinguish the genuine DNA fragments from those of modern contaminants and from those of microorganisms that are co-extracted from the sample
Since ancient DNA has characteristic damage patterns, e.g. a high degree of cytosine deaminations on the five prime ends, and at the three prime end, there is a high frequency of Gs changing to As.
Match these ancient DNA fragments to a map of the complete human genome reference sequence; this allowed the book of the target's genome to be assembled from the millions of random sequence fragments recovered from the bones
The more times we had a read confirmed by multiple fragments, the more confident the sequencing could be in distinguishing the true DNA sequence, it is the sequencing depth issue
Then genotype samples to discover, filter, and analyze variants.


text:
How to convert VCF files to input file format required by XP-CLR1.0 software?
I attempt to do a selective sweep analysis by using the software XP-CLR1.0. My genotypic data is VCF format, so I wonder if someone with experience on this package could recommend any tool to convert from VCF to the input file format required by this software?
Output files I need are:

geno file: contains the genotype of a single SNP, for example:

1 0 1 1 9 9
1 1 1 0 0 0

It contains two SNPs (rows) from three individuals (columns). First two columns are the two alleles of individual #1. Values 1 and 0 correspond to two allele types. 9 denotes missing data.

snp file: each row contains information of each SNP, for example:

snp_name, chr#, genetic_distance(cM), Physical_distance, Ref_allele, Alt_allele

Please [edit] your question and explain exactly what the expected output format is. Ideally, add a small example of your input file and then the output you would expect from that example. Read the tool's documentation, read the paper describing it, understand the format and then tell us what it is and we can help you get there.

Please let me know if the question clear now.

Not really, I'm afraid. What is the expected output from your example input? What is the exact input and output format you have and need?

What I need to do is to convert VCF files to both of the formats explained above in my question (geno.file and snp.file).

OK, so please add an example input, one or two variant lines from a VCF file, and then show the output you would want from that input. Tell us how we can know what values "genetic_distance", "physical_distance" etc., should have, where that information is in the vcf file.

Perhaps https://github.com/hardingnj/xpclr would solve the problem? Looks like this version of the software runs the XP-CLR algorithm and accepts vcf files as input. Probably worth a look


text:
How to calculate Allele frequency from vcf file
Organism under investigation is  Plasmodium falciparum. How to calculate the allele frequency for each row?

I tried with this code:
    import pandas as pd

# Read CSV file
# df = pd.read_csv("Pf_M76611.pf7.200_samples - Copy.csv")

df= pd.read_csv("/content/Pf_M76611.pf7.200_samples - Copy.csv")

# Define function to calculate allele frequencies
def calc_allele_freq(row):
    # Get total number of samples
    num_samples = len(row) - 9

    # Get reference allele
    ref_allele = row['REF']

    # Get alternative alleles
    alt_alleles = row['ALT'].split(',')

    # Initialize list to store allele frequencies
    allele_freqs = []

    # Loop over alternative alleles
    for alt_allele in alt_alleles:
        # Initialize count for this allele
        count = 0

        # Loop over samples
        for i in range(9, len(row)):
            # Get genotype for this sample
            gt = row[i].split(':')[0]

            # Count this allele if it is present in the genotype
            if alt_allele in gt:
                count += 1

        # Calculate frequency for this allele
        freq = count / num_samples

        # Add frequency to list
        allele_freqs.append(freq)

    # Add frequency for reference allele
    allele_freqs.insert(0, 1 - sum(allele_freqs))

    # Return allele list and frequency list as a tuple
    allele_list = [ref_allele] + alt_alleles
    freq_list = [round(freq, 1) for freq in allele_freqs]
    return (allele_list, freq_list)

# Apply function to each row of the dataframe
allele_freq_df = df.apply(calc_allele_freq, axis=1, result_type='expand')

# Add columns to dataframe
df['Allele List'] = allele_freq_df[0]
df['Proportion List'] = allele_freq_df[1]

# Select columns for output and drop duplicates
# output_df = df[['#CHROM', 'POS', 'Allele List', 'Proportion List']].drop_duplicates()

# # Convert output dataframe to string and print
# print(output_df.to_string(index=False))

Could you please show that table as text format? It's quite difficult to read on my mobile device.

Coding: first and foremost all operations should be performed directly in pandas. It is much more efficient and alot easier. You don't need all the loops, pandas does that. apply is okay. I would use a very different strategy (below).

The actually contents of the csv are not clear.
This is what I think might have happened: there's   multiple VCFs against the same reference I assume these were loaded and a database join operation performed using the reference genome column 'REF' as the index for the join across all input VCFs. I've placed this here because I'm not 100% sure what the input file represents - if it's not what I think it is - thats how I'd get to the allelic frequency.
Why do I think that?

Without multiple strains the allelic frequency cannot be calculated.
The file states '200 samples'
The presence of multiple alleles in the 'ALT' column.

The problems are:

heterozygosity. This assumes falciparum is in its haploid state (it has a diploid stage). If its the diploid stage the allelic frequency cannot be calculated using this input file because the linkage was discarded (which is really important)
The underlying problem is the number of actual strains isn't explicitly stated .
the frequency of the 'ALT' doesn't appear to be declared
even if the a frequency was declared (but not presented) that must be linked to the 'ALT' column, each allele needs a exact frequency, otherwise population genetics will not work on this data.

First conclusion The allelic frequency is performed across each row using the reference genome as the reference via value_counts() divided by the number of strains and presented in its own column. Somehow the relevant columns are missing - will need to see them.
However Again on the information presented I don't think it can be done because if this represents 200 samples the frequency of a given allele is not declared in this data, or simply the number of 'ALT's across all VCFs. Without accessing that number the allelic frequency is not doable because it's entirely dependent on it.
Possible alternative conclusion Its possible the question needs to be redefined not as allelic frequency - this has an exact definition within population genetics, but genetic diversity, i.e. remove the term 'frequency' altogether. What I think is being calculated is genetic diversity of all samples in the file against a reference genome. If that is correct it would be good to change the question accordingly. These are different mathematical entities.
Whatever operation is required is easy to calculate, but its obtaining that information from the input file and the precise nature of the investigation.


text:
Singularity & GATK 4.4.0.0. CreateSequenceDictionary
I ran the following code in bash (I added print statements for ease):
# Module loading

echo "I am loading the specified modules"

module load Singularity/version_number

# Define the path to the singularity:

singularity_image="/path/to/folder/GATK/gatk_latest.sif"

echo "I am going to make the `.dict ` file"

singularity exec $singularity_image /gatk/gatk --java-options -Xms8000m \
  CreateSequenceDictionary \
  -R /path/to/my/primary_assembly.fa \
  -O /path/to/the/primary_assembly.dict

echo "I am all done"

However, I get the following out file:
I am loading the specified modules
job_scripts_number: line 23:
.dict: command not found I am going to make the  file
INFO  NativeLibraryLoader - Loading libgkl_compression.so from
jar:file:/gatk/gatk-package-4.4.0.0-local.jar!/path/to/libgkl_compression.so [Today's date] CreateSequenceDictionary --OUTPUT
path/to/my/primary_assembly.dict
--REFERENCE /path/to/my/primary_assembly.fa
--TRUNCATE_NAMES_AT_WHITESPACE true --NUM_SEQUENCES 2147483647 --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 999999 --CREATE_INDEX false --CREATE_MD5_FILE false --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false [Today's date] Executing as user@node_number on Linux version
; OpenJDK 64-Bit Server VM version;
Deflater: A; Inflater: A; Provider is available; Picard
version: Version:4.4.0.0 [Today's date]
picard.sam.CreateSequenceDictionary done. Elapsed time: 0.53 minutes.
Runtime.totalMemory()=8388608000 Tool returned: 0 Using GATK jar
/gatk/gatk-package-4.4.0.0-local.jar Running:
    java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -Xms8000m -jar /gatk/gatk-package-4.4.0.0-local.jar CreateSequenceDictionary -R
/path/to/primary_assembly.fa
-O /path/to/primary_assembly.dict
I am all done

Why would I get the error ".dict: command not found I am going to make the  file " especially when I look at the head of the .dict file and it seems to look the way it should?

To make debugging easier, ensure you include #!/bin/bash and set -eo pipefail as the first and second lines of your script. The former specifies the interpreter and the latter fails the entire script the moment one command fails. If you'd had that line, your script would have stopped executing at .dict: command not found and the last line of your STDERR would tell you where the problem was.

This has nothing to do with Singularity/Docker or GATK. You get the error .dict: command not found because you've used backticks to perform a command substitution in a prior echo command:

echo "I am going to make the `.dict ` file"

Backticks are legacy syntax, you could have also written:
echo "I am going to make the $(.dict) file"

But neither of these makes any sense - there is no .dict command. Of course single quotes would prevent the command substitution:
echo 'I am going to make the `.dict` file'

But this might cause some confusion when the output is displayed. A better choice would be to use double quotes, but without the backticks. For example:
echo "I am going to make the '.dict' file"

It worked! Thank you I used the echo "I am going to make the '.dict' file"