text:
How can I retrieve the GRCh38 coordinates of a list of rsids?
I have a list of about 100,000 rsids and I want to get their genomic coordinates on the GRCh38 genome build. Is there a command line tool that allows me to do this? If yes, which one?
I have tried bcftools but, given the error message I got, I believe that it does not support search by rsid:
bcftools view -R rsids_chr1_subset.txt -f PASS -v snps dbSNP_v156_GRCh38.p14.vcf.gz > rsids_chr1.vcf
[E::bcf_sr_regions_init] Could not parse 1-th line of file rsids_chr1_subset.txt, using the columns 1,2[,-1]
Failed to read the regions: rsids_chr1_subset.txt

The tricky thing is doing this fast. If I remember correctly, I have successfully used an implementation written in C a long time ago that did just this. It also had a Python version. The name escapes me at the moment. However, this just reminded me that I ended up writing my own implementation because there was an issue. I will make an answer.

Several years ago I created rsidx to address this particular use case. https://github.com/bioforensics/rsidx
# VCF should be sorted by genomic coordinates and indexed by tabix
rsidx index dbSNP151_GRCh38.vcf.gz dbSNP151_GRCh38.rsidx
rsidx search dbSNP151_GRCh38.vcf.gz dbSNP151_GRCh38.rsidx rs3114908 rs10756819

rsids_chr1.txt contains a list of rsids, one per line.
bcftools view --include ID==@rsids_chr1.txt dbSNP_v156_GRCh38.p14.vcf.gz -Ou | bcftools query -f '%CHROM  %ID %POS  %REF  %ALT{0}\n' > dbSNP_v156_GRCh38.p14_AF.vcf

This was part of a custom liftover function which handled multiple alleles and reverse compliments. I do not know if it will be useful. I believe there was a separate thing somewhere with lookup that I may add. It handled apparent non-1:1 mappings in addition to platform-specific "rsIDs."
converter = get_lifter(current, target)
def convert_coordinates(input_file, output_file, converter):
        line_count = 0
        skipped = 0
        strand_flips = 0
        with open(input_file, 'r') as infile, open(output_file, 'w') as outfile:
            for line in infile:
                # Skip headers
                if line.startswith("#"):
                    outfile.write(line)
                    continue

                # VCF files are tab-separated
                parts = line.strip().split('\t')

                # Get the chromosome and position
                chrom = parts[0]
                pos = int(parts[1])

                # Convert the coordinates
                try:
                    new_coords = converter.query(chrom, pos)
                except KeyError:
                    #print(f"Chromosome {chrom} not found in the liftover chain file.")
                    skipped += 1
                    continue

                # If there's no new coordinates, we skip this variant
                if not new_coords:
                    skipped += 1
                    continue

                # Get the new chromosome, position and strand
                new_chrom, new_pos, strand = new_coords[0]

                # If the new strand is negative, reverse complement the REF and ALT alleles
                if strand == '-':
                    parts[3] = reverse_complement(parts[3])
                    parts[4] = reverse_complement(parts[4])
                    strand_flips += 1

                # Replace the old chromosome and position with the new ones
                parts[0] = new_chrom
                parts[1] = str(new_pos)

                # Write the line back into the file
                outfile.write('\t'.join(parts))
                outfile.write('\n')

                # Print a progress statement every 1000 lines
                line_count += 1
                if line_count % 100000 == 0:
                    print(f"Processed {line_count} lines")

        # Calculate and print the percentage of variants skipped
        percentage_skipped = (skipped / line_count) * 100
        percentage_flips = (strand_flips / line_count) * 100
        print(f"\nFinished processing. Out of {line_count} total lines:")
        print(f"{percentage_skipped}% of variants were skipped due to missing chromosomal data.")
        print(f"{percentage_flips}% of variants required a strand flip.")


text:
which secondary structure assignment algorithm is used by PDB?
older releases of PDB used to specify DSSP as the 2ary structure assignment algorithm used in PDB, right now no algorithm is explicitly stated. Where can I find a reference for it in the database?

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.

Does the PDB assign the secondary structure tags, or just validate the one uploaded in the repository ??

DSSP 4.0 ?? https://github.com/PDB-REDO/dssp is used in PDB-redo

cant find any reference tough

source: Stack Exchange

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

sort: votes


text:
How to generate rsid values locally from VCF file
I recently used Dante to obtain a WGS and have downloaded several files produced during this process. There is a file that ends '.vcf' which I did think was a VCF file but I'm not sure now after a previous question. Part of that file is included below. Its about 1.1GB in size.
It appears to me to follow VCF format at least to a degree but it is missing the rsid's. The field where an rsid would usually be in a VCF file is replaced with a period (.).
I would like to backfill (annotate?) this file with known rsid's (where defined). There are nearly 4 million of them and online tools don't seem ideal, I would prefer to be able to do it locally, if that is possible. I am not looking for software to do this, but a data file which specifies which rsid is assigned to a particular chr:pos:allele. That may or may not be possible - perhaps I misunderstand and its far more complex than that. But from what I can tell, for a particular reference genome, and a particular build of dbsnp, there is a simple relationship as described.
I want to do this, for instance, to be able to use various tools, like SNPedia. Dante appears to use GRCh37, vs SNPedia using mainly GRCh38, which I can see is an issue, and will probably require me to buy a GRCh38 conversion from Dante (or do it myself). But please ignore that particular issue right now.
I don't quite understand what the dependencies are. I've seen comments that suggest that reference genome is not important to generating rsid's from VCF files and that only dbsnp build number needs to match the environment that the VCF file was generated in, but that does not seem to make much sense to me. From what I can see, if you are converting chr:pos information to rsid's, you need to be using the right reference genome. Does dbsnp build really matter? Should you not use the latest build, if it is continually being added to - surely rsid's are not reused? And then a post on the ANNOVAR site seems to suggest that no one really knows what rsid's are and that to different people they are different things, further muddying the waters.
So in terms of question - is there a file, downloadable from somewhere, that contains mappings from chr:pos to rsid, and if so, how does one select the correct file for a particular dataset (IE generated with a particular reference genome, with a particular dbsnp, to be used with references that have used a particular reference/dbsnp build like GRCh38 / dbsnp build 151)?
This file (header first), downloaded from Dante, is named .filtered.snp.vcf
    ##fileformat=VCFv4.2
    ##FILTER=<ID=PASS,Description="All filters passed">
    ##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths (counting only informative reads out of the total reads) for the ref and alt alleles in the order listed">
    ##FORMAT=<ID=AF,Number=A,Type=Float,Description="Allele fractions for alt alleles in the order listed">
    ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
    ##FORMAT=<ID=F1R2,Number=R,Type=Integer,Description="Count of reads in F1R2 pair orientation supporting each allele">
    ##FORMAT=<ID=F2R1,Number=R,Type=Integer,Description="Count of reads in F2R1 pair orientation supporting each allele">
    ##FORMAT=<ID=GP,Number=G,Type=Float,Description="Phred-scaled posterior probabilities for genotypes as defined in the VCF specification">
    ##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
    ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
    ##FORMAT=<ID=MB,Number=4,Type=Integer,Description="Per-sample component statistics to detect mate bias">
    ##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">
    ##FORMAT=<ID=PRI,Number=G,Type=Float,Description="Phred-scaled prior probabilities for genotypes">
    ##FORMAT=<ID=PS,Number=1,Type=Integer,Description="Physical phasing ID information, where each unique ID within a given sample (but not across samples) connects records within a phasing group">
    ##FORMAT=<ID=SB,Number=4,Type=Integer,Description="Per-sample component statistics which comprise the Fisher's Exact Test to detect strand bias">
    ##FORMAT=<ID=SQ,Number=A,Type=Float,Description="Somatic quality">
    ##DRAGENCommandLine=<ID=HashTableBuild,Version="SW: 01.003.044.3.9.5, HashTableVersion: 8",CommandLineOptions="dragen --build-hash-table true --output-directory grch37 --ht-reference grch37.fasta">
    ##DRAGENCommandLine=<ID=dragen,Version="SW: 05.121.645.4.0.3, HW: 05.121.645",Date="Thu Mar 02 09:58:00 UTC 2023",CommandLineOptions="-f -c output.cfg -r /references/grch37">
    ##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">
    ##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed">
    ##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
    ##INFO=<ID=FS,Number=1,Type=Float,Description="Phred-scaled p-value using Fisher's exact test to detect strand bias">
    ##INFO=<ID=QD,Number=1,Type=Float,Description="Variant Confidence/Quality by Depth">
    ##INFO=<ID=SOR,Number=1,Type=Float,Description="Symmetric Odds Ratio of 2x2 contingency table to detect strand bias">
    ##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (informative and non-informative); some reads may have been filtered based on mapq etc.">
    ##INFO=<ID=FractionInformativeReads,Number=1,Type=Float,Description="The fraction of informative reads out of the total reads">
    ##INFO=<ID=MQ,Number=1,Type=Float,Description="RMS Mapping Quality">
    ##INFO=<ID=MQRankSum,Number=1,Type=Float,Description="Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities">
    ##INFO=<ID=ReadPosRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias">
    ##FILTER=<ID=DRAGENSnpHardQUAL,Description="Set if true:QUAL < 3">
    ##FILTER=<ID=DRAGENIndelHardQUAL,Description="Set if true:QUAL < 3">
    ##FILTER=<ID=LowDepth,Description="Set if true:DP <= 1">
    ##FILTER=<ID=PloidyConflict,Description="Genotype call from variant caller not consistent with chromosome ploidy">
    ##FILTER=<ID=base_quality,Description="Site filtered because median base quality of alt reads at this locus does not meet threshold">
    ##FILTER=<ID=filtered_reads,Description="Site filtered because too large a fraction of reads have been filtered out">
    ##FILTER=<ID=fragment_length,Description="Site filtered because absolute difference between the median fragment length of alt reads and median fragment length of ref reads at this locus exceeds threshold">
    ##FILTER=<ID=low_af,Description="Allele frequency does not meet threshold">
    ##FILTER=<ID=low_depth,Description="Site filtered because the read depth is too low">
    ##FILTER=<ID=low_frac_info_reads,Description="Site filtered because the fraction of informative reads is below threshold">
    ##FILTER=<ID=low_normal_depth,Description="Site filtered because the normal sample read depth is too low">
    ##FILTER=<ID=long_indel,Description="Site filtered because the indel length is too long">
    ##FILTER=<ID=mapping_quality,Description="Site filtered because median mapping quality of alt reads at this locus does not meet threshold">
    ##FILTER=<ID=multiallelic,Description="Site filtered because more than two alt alleles pass tumor LOD">
    ##FILTER=<ID=non_homref_normal,Description="Site filtered because the normal sample genotype is not homozygous reference">
    ##FILTER=<ID=no_reliable_supporting_read,Description="Site filtered because no reliable supporting somatic read exists">
    ##FILTER=<ID=panel_of_normals,Description="Seen in at least one sample in the panel of normals vcf">
    ##FILTER=<ID=read_position,Description="Site filtered because median of distances between start/end of read and this locus is below threshold">
    ##FILTER=<ID=RMxNRepeatRegion,Description="Site filtered because all or part of the variant allele is a repeat of the reference">
    ##FILTER=<ID=str_contraction,Description="Site filtered due to suspected PCR error where the alt allele is one repeat unit less than the reference">
    ##FILTER=<ID=too_few_supporting_reads,Description="Site filtered because there are too few supporting reads in the tumor sample">
    ##FILTER=<ID=weak_evidence,Description="Somatic variant score does not meet threshold">
    ##ML=PackageVersion=4,ModelVersion=4,ModelCRC=2921457924
    ##contig=<ID=1,length=249250621>
    ##contig=<ID=2,length=243199373>
    ##contig=<ID=3,length=198022430>
    ##contig=<ID=4,length=191154276>
    ##contig=<ID=5,length=180915260>
    ##contig=<ID=6,length=171115067>
    ##contig=<ID=7,length=159138663>
    ##contig=<ID=8,length=146364022>
    ##contig=<ID=9,length=141213431>
    ##contig=<ID=10,length=135534747>
    ##contig=<ID=11,length=135006516>
    ##contig=<ID=12,length=133851895>
    ##contig=<ID=13,length=115169878>
    ##contig=<ID=14,length=107349540>
    ##contig=<ID=15,length=102531392>
    ##contig=<ID=16,length=90354753>
    ##contig=<ID=17,length=81195210>
    ##contig=<ID=18,length=78077248>
    ##contig=<ID=19,length=59128983>
    ##contig=<ID=20,length=63025520>
    ##contig=<ID=21,length=48129895>
    ##contig=<ID=22,length=51304566>
    ##contig=<ID=X,length=155270560>
    ##contig=<ID=Y,length=59373566>
    ##contig=<ID=MT,length=16569>
    ##contig=<ID=GL000207.1,length=4262>
    ##contig=<ID=GL000226.1,length=15008>
    ##contig=<ID=GL000229.1,length=19913>
    ##contig=<ID=GL000231.1,length=27386>
    ##contig=<ID=GL000210.1,length=27682>
    ##contig=<ID=GL000239.1,length=33824>
    ##contig=<ID=GL000235.1,length=34474>
    ##contig=<ID=GL000201.1,length=36148>
    ##contig=<ID=GL000247.1,length=36422>
    ##contig=<ID=GL000245.1,length=36651>
    ##contig=<ID=GL000197.1,length=37175>
    ##contig=<ID=GL000203.1,length=37498>
    ##contig=<ID=GL000246.1,length=38154>
    ##contig=<ID=GL000249.1,length=38502>
    ##contig=<ID=GL000196.1,length=38914>
    ##contig=<ID=GL000248.1,length=39786>
    ##contig=<ID=GL000244.1,length=39929>
    ##contig=<ID=GL000238.1,length=39939>
    ##contig=<ID=GL000202.1,length=40103>
    ##contig=<ID=GL000234.1,length=40531>
    ##contig=<ID=GL000232.1,length=40652>
    ##contig=<ID=GL000206.1,length=41001>
    ##contig=<ID=GL000240.1,length=41933>
    ##contig=<ID=GL000236.1,length=41934>
    ##contig=<ID=GL000241.1,length=42152>
    ##contig=<ID=GL000243.1,length=43341>
    ##contig=<ID=GL000242.1,length=43523>
    ##contig=<ID=GL000230.1,length=43691>
    ##contig=<ID=GL000237.1,length=45867>
    ##contig=<ID=GL000233.1,length=45941>
    ##contig=<ID=GL000204.1,length=81310>
    ##contig=<ID=GL000198.1,length=90085>
    ##contig=<ID=GL000208.1,length=92689>
    ##contig=<ID=GL000191.1,length=106433>
    ##contig=<ID=GL000227.1,length=128374>
    ##contig=<ID=GL000228.1,length=129120>
    ##contig=<ID=GL000214.1,length=137718>
    ##contig=<ID=GL000221.1,length=155397>
    ##contig=<ID=GL000209.1,length=159169>
    ##contig=<ID=GL000218.1,length=161147>
    ##contig=<ID=GL000220.1,length=161802>
    ##contig=<ID=GL000213.1,length=164239>
    ##contig=<ID=GL000211.1,length=166566>
    ##contig=<ID=GL000199.1,length=169874>
    ##contig=<ID=GL000217.1,length=172149>
    ##contig=<ID=GL000216.1,length=172294>
    ##contig=<ID=GL000215.1,length=172545>
    ##contig=<ID=GL000205.1,length=174588>
    ##contig=<ID=GL000219.1,length=179198>
    ##contig=<ID=GL000224.1,length=179693>
    ##contig=<ID=GL000223.1,length=180455>
    ##contig=<ID=GL000195.1,length=182896>
    ##contig=<ID=GL000212.1,length=186858>
    ##contig=<ID=GL000222.1,length=186861>
    ##contig=<ID=GL000200.1,length=187035>
    ##contig=<ID=GL000193.1,length=189789>
    ##contig=<ID=GL000194.1,length=191469>
    ##contig=<ID=GL000225.1,length=211173>
    ##contig=<ID=GL000192.1,length=547496>
    ##contig=<ID=NC_007605,length=171823>
    ##contig=<ID=hs37d5,length=35477943>
    ##reference=file:///references/grch37/reference.bin
    ##bcftools_filterVersion=1.9+htslib-1.9
    ##bcftools_filterCommand=filter -i 'TYPE="snp" ' /processing/sample.raw.vcf.gz; Date=Thu Mar  2 10:28:45 2023
    #CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  GFX0xxxxxx
    1   12783   .   G   A   6.16    PASS    AC=1;AF=0.5;AN=2;DP=37;FS=0;MQ=11.6;MQRankSum=-2.137;QD=0.45;ReadPosRankSum=2.878;SOR=1.949;FractionInformativeReads=0.973  GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:10,26:0.722:36:4,14:6,12:4:38,0,1:6.1614,2.5231,7.0198:0,34.77,37.77:10,0,26,0:4,6,14,12
    1   13868   .   A   G   4.69    PASS    AC=1;AF=0.5;AN=2;DP=61;FS=2.275;MQ=16.39;MQRankSum=2.368;QD=0.36;ReadPosRankSum=-1.078;SOR=1.136;FractionInformativeReads=1 GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:36,25:0.41:61:16,15:20,10:4:37,0,8:4.6918,2.1112,13.42:0,34.77,37.78:15,21,8,17:23,13,12,13
    1   13980   .   T   C   23.14   PASS    AC=2;AF=1;AN=2;DP=2;FS=0;MQ=5.7;MQRankSum=0;QD=6.5;ReadPosRankSum=0;SOR=2.303;FractionInformativeReads=1    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    1/1:0,2:1:2:0,1:0,1:19:61,24,0:23.143,20.887,0.056833:0,34.77,37.77:0,0,2,0:0,0,1,1
    1   14464   .   A   T   14.11   PASS    AC=1;AF=0.5;AN=2;DP=48;FS=7.243;MQ=22.38;MQRankSum=3.14;QD=0.86;ReadPosRankSum=2.527;SOR=0.267;FractionInformativeReads=1   GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:11,37:0.771:48:8,16:3,21:7:77,6,0:14.113,0.98271,7.8591:0,34.77,37.77:3,8,19,18:4,7,20,17
    1   14907   .   A   G   3.67    PASS    AC=1;AF=0.5;AN=2;DP=77;FS=0.893;MQ=26.78;MQRankSum=-0.058;QD=0.41;ReadPosRankSum=5.956;SOR=0.532;FractionInformativeReads=1 GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:48,29:0.377:77:31,13:17,16:4:36,0,22:3.6675,2.453,27.555:0,34.77,37.77:27,21,15,14:28,20,15,14
    1   14930   .   A   G   3.26    PASS    AC=1;AF=0.5;AN=2;DP=82;FS=3.051;MQ=27.04;MQRankSum=-0.363;QD=0.37;ReadPosRankSum=-0.277;SOR=0.754;FractionInformativeReads=1    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:51,31:0.378:82:30,16:21,15:3:35,0,21:3.2627,2.7897,26.606:0,34.77,37.77:28,23,14,17:29,22,17,14
    1   15211   .   T   G   6.66    PASS    AC=1;AF=0.5;AN=2;DP=70;FS=0;MQ=21.21;MQRankSum=0.533;QD=0.49;ReadPosRankSum=-0.905;SOR=0.743;FractionInformativeReads=1 GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:28,42:0.6:70:12,22:16,20:6:40,0,11:6.6596,1.2235,15.27:0,34.77,37.77:15,13,23,19:15,13,17,25

Then a couple of chunks from random bits of the file:
offset 100000 lines:
    1   82279586    .   T   C   62.12   PASS    AC=2;AF=1;AN=2;DP=25;FS=0;MQ=250;MQRankSum=0;QD=4.89;ReadPosRankSum=0;SOR=0.941;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB  1/1:0,25:1:25:0,10:0,15:52:100,55,0:62.124,52.173,2.8992e-05:0,34.77,37.77:0,0,11,14:0,0,10,15
    1   82280153    .   G   A   62.73   PASS    AC=2;AF=1;AN=2;DP=31;FS=0;MQ=250;MQRankSum=0;QD=4.53;ReadPosRankSum=0;SOR=1.044;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB  1/1:0,31:1:31:0,15:0,16:57:101,61,0:62.731,57.663,9.319e-06:0,34.77,37.77:0,0,18,13:0,0,14,17
    1   82280683    .   C   T   54.83   PASS    AC=2;AF=1;AN=2;DP=17;FS=0;MQ=250;MQRankSum=0;QD=5.77;ReadPosRankSum=0;SOR=1.061;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB  1/1:0,17:1:17:0,8:0,9:45:93,49,0:54.832,45.914,0.00012581:0,34.77,37.77:0,0,10,7:0,0,11,6
    1   82280870    .   T   C   41.24   PASS    AC=1;AF=0.5;AN=2;DP=19;FS=1.808;MQ=250;MQRankSum=3.592;QD=2.58;ReadPosRankSum=1.858;SOR=1.14;FractionInformativeReads=1 GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:11,8:0.421:19:8,2:3,6:40:76,0,45:41.236,0.00039866,47.817:0,34.77,37.77:4,7,2,6:5,6,6,2
    1   82280975    .   A   T   40.52   PASS    AC=1;AF=0.5;AN=2;DP=24;FS=1.569;MQ=250;MQRankSum=4.07;QD=2.02;ReadPosRankSum=1.025;SOR=1.085;FractionInformativeReads=1 GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:14,10:0.417:24:5,8:9,2:40:75,0,50:40.524,0.00040643,53.142:0,34.77,37.77:7,7,4,6:8,6,4,6
    1   82281492    .   T   C   44.14   PASS    AC=1;AF=0.5;AN=2;DP=22;FS=6.923;MQ=250;MQRankSum=3.873;QD=2.2;ReadPosRankSum=-0.334;SOR=2.064;FractionInformativeReads=1    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:13,9:0.409:22:9,4:4,5:44:79,0,51:44.143,0.00018328,54.3:0,34.77,37.77:6,7,7,2:6,7,6,3
    1   82282062    .   T   C   8.34    PASS    AC=1;AF=0.5;AN=2;DP=26;FS=1.974;MQ=247.84;MQRankSum=3.383;QD=1.04;ReadPosRankSum=-1.561;SOR=0.863;FractionInformativeReads=1    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:21,5:0.192:26:9,3:12,2:8:42,0,29:8.3401,0.69098,32.665:0,34.77,37.77:9,12,3,2:9,12,1,4
    1   82282378    .   C   T   42.61   PASS    AC=1;AF=0.5;AN=2;DP=24;FS=1.685;MQ=250;MQRankSum=4.07;QD=2.02;ReadPosRankSum=2.898;SOR=1.23;FractionInformativeReads=1  GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:14,10:0.417:24:8,6:6,4:43:77,0,57:42.606,0.0002423,60.191:0,34.77,37.77:8,6,7,3:10,4,7,3
    1   82282788    .   C   T   48.99   PASS    AC=1;AF=0.5;AN=2;DP=30;FS=0;MQ=245.71;MQRankSum=2.889;QD=1.67;ReadPosRankSum=1.975;SOR=0.914;FractionInformativeReads=1 GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:14,16:0.533:30:6,9:8,7:46:84,0,46:48.985,0.00010665,49.217:0,34.77,37.77:8,6,10,6:6,8,7,9
    1   82282947    .   G   T   47.31   PASS    AC=1;AF=0.5;AN=2;DP=32;FS=10.606;MQ=250;MQRankSum=4.796;QD=1.56;ReadPosRankSum=0.68;SOR=0.725;FractionInformativeReads=1    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:15,17:0.531:32:6,8:9,9:46:82,0,48:47.312,0.00011286,51.287:0,34.77,37.77:11,4,7,10:8,7,8,9

offset 2000000 lines:
    8   142260995   .   T   C   65.76   PASS    AC=2;AF=1;AN=2;DP=28;FS=0;MQ=248.93;MQRankSum=0;QD=4.69;ReadPosRankSum=0;SOR=1.179;FractionInformativeReads=1   GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    1/1:0,28:1:28:0,14:0,14:53:104,56,0:65.759,52.946,2.3297e-05:0,34.77,37.77:0,0,11,17:0,0,13,15
    8   142261389   .   T   C   66.03   PASS    AC=2;AF=1;AN=2;DP=28;FS=0;MQ=250;MQRankSum=0;QD=4.67;ReadPosRankSum=0;SOR=0.997;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB  1/1:0,28:1:28:0,12:0,16:54:104,57,0:66.026,54.004,1.812e-05:0,34.77,37.77:0,0,16,12:0,0,14,14
    8   142261541   .   G   A   65.34   PASS    AC=2;AF=1;AN=2;DP=30;FS=0;MQ=250;MQRankSum=0;QD=4.51;ReadPosRankSum=0;SOR=1.765;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB  1/1:0,30:1:30:0,14:0,16:52:103,55,0:65.339,52.429,2.5886e-05:0,34.77,37.77:0,0,9,21:0,0,17,13
    8   142261618   .   A   G   64.46   PASS    AC=2;AF=1;AN=2;DP=27;FS=0;MQ=250;MQRankSum=0;QD=4.57;ReadPosRankSum=0;SOR=1.302;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB:PS   1|1:0,27:1:27:0,14:0,13:53:102,56,0:64.465,53.434,2.1227e-05:0,34.77,37.78:0,0,10,17:0,0,15,12:142261618
    8   142261621   .   C   T   60.33   PASS    AC=2;AF=1;AN=2;DP=25;FS=0;MQ=250;MQRankSum=0;QD=4.7;ReadPosRankSum=0;SOR=1.051;FractionInformativeReads=0.96    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB:PS 1|1:0,24:1:24:0,13:0,11:48:98,51,0:60.326,48.044,7.2481e-05:0,34.77,37.78:0,0,10,14:0,0,13,11:142261618
    8   142261967   .   C   T   68.75   PASS    AC=2;AF=1;AN=2;DP=37;FS=0;MQ=250;MQRankSum=0;QD=4.26;ReadPosRankSum=0;SOR=0.859;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB  1/1:0,37:1:37:0,17:0,20:56:107,60,0:68.745,56.755,9.8367e-06:0,34.77,37.77:0,0,20,17:0,0,18,19
    8   142262327   .   A   G   68.48   PASS    AC=2;AF=1;AN=2;DP=42;FS=0;MQ=250;MQRankSum=0;QD=4.13;ReadPosRankSum=0;SOR=0.788;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB  1/1:0,42:1:42:0,21:0,21:56:106,60,0:68.476,56.786,9.8367e-06:0,34.77,37.77:0,0,22,20:0,0,11,31
    8   142262413   .   G   A   62.88   PASS    AC=2;AF=1;AN=2;DP=27;FS=0;MQ=231.86;MQRankSum=0;QD=4.75;ReadPosRankSum=0;SOR=1.302;FractionInformativeReads=1   GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    1/1:0,27:1:27:0,16:0,11:51:101,54,0:62.881,51.183,3.5205e-05:0,34.77,37.77:0,0,10,17:0,0,9,18
    8   142263074   .   G   A   51.24   PASS    AC=2;AF=1;AN=2;DP=20;FS=0;MQ=143.27;MQRankSum=0;QD=5.31;ReadPosRankSum=0;SOR=0.892;FractionInformativeReads=1   GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    1/1:0,20:1:20:0,9:0,11:36:89,39,0:51.237,35.702,0.0012015:0,34.77,37.77:0,0,11,9:0,0,10,10
    8   142263171   .   C   T   52.9    PASS    AC=2;AF=1;AN=2;DP=22;FS=0;MQ=122.2;MQRankSum=0;QD=4.81;ReadPosRankSum=0;SOR=0.693;FractionInformativeReads=0.909    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    1/1:0,20:1:20:0,11:0,9:35:91,38,0:52.899,34.894,0.0014297:0,34.77,37.77:0,0,10,10:0,0,9,11

Download the dbSNP VCF file (latest version) and map to your input VCF using bcftools annotate. You'll want to match reference sequence version (hg19/hg38) and use the latest version of dbSNP to get the most up to date information.
dbSNP VCF FTP folder: https://ftp.ncbi.nih.gov/snp/latest_release/VCF/

Thanks, I eventually got this working although the brittleness with which it treats input file compression standards is rather annoying and took me ages to work out, hence my delay. I've finally annotated the vcf with rsid with bcftools, so thanks very much. Since my data is refed to GRCh37, and SNPedia is made with 38, although it works I wonder if it's worth 'getting' 38 somehow. Looks like generating it from the BAM or fastq would be the cleanest way to do this. Could I ask you to comment on that, and recommend any software or guides?

Yes, re-aligning to hg38 is the way to go. liftOver is an option too, but it has its shortcomings.

So regarding the realigning, is that a massive undertaking or relatively simple? I imagine it's simple - just testing the original files (bam or fastq etc) against the new reference genome hg38. Or does it require detailed info about the original sequencing operation that I am not likely to be able to get out of Dante (as in gringers reply re illumina snpchip versions being matched with specific dbsnp builds etc)?

Since it's WGS, my guess is that realigning should work just fine. Targeted sequencing might have created problems. You can ask @gringer for their take.

RamRS - thanks. Gringer, wonder if you come by these parts again if you could comment? As an extra ask, hoping its not too close to the wind to ask on SE about software recommendations, any suggestions for a tool to do this (I'm on Linux so Windows only software is out for me)?

As a partial answer to the specific questions about definitions:

rsIDs define a specific location and type of variation*, but the location of that variation can change within different reference genomes. This may include a change in chromosome, as well as a change in position.

Different dbSNP builds can add additional variation, but the same rsID will always refer to the same variant in all dbSNP builds that refer to that variation.

rsIDs can be invalidated in new dbSNP builds for various reasons. It's a good idea to always use the latest dbSNP build when annotating variants.

rsIDs are never recycled. There is a prescribed process for merging and splitting variants that ensures backwards compatibility. If a variation is merged, I think the approach is to always define the earliest-numbered variant as the primary variation. If a variation is split, I think the approach is to use a new rsID for each of the separate variations.

* by "variation", I'm referring to a type of modification in the DNA sequence, rather than a specific instance of that modification. This is explained in the dbSNP documentation (an "ss" is a submitted variant, a submission made to dbSNP):

Each ss that maps to the same position and is the same variant type is assigned to an existing RefSNP (rs) or assigned a new one.

Although that isn't entirely clear, this suggests that if a new variant is submitted that falls at the same position of an exiting one and is of the same type as the existing one, then the new one may simply be added to the same rsID as the existing one. So, for example, if we have a single nucleotide G => C variant at position N, and a new G => T variant is reported at the same position, then the new one could have the same rsID as the old one.

In response to an additional comment / question:

a post on the ANNOVAR site seems to suggest that no one really knows what rsid's are and that to different people they are different things, further muddying the waters.... Is this guy an authority on this or not?

I can't comment on authority for that guy; the authorship is unclear to me. What matters is how the dbSNP people understand dbSNP (i.e. the documentation link above), because they are the authority on it.
My understanding of it was helped because many years ago I attended a close-knit workshop that included a talk by Mike Feolo, and I had been working on a human dataset with many oddities, so I was able to ask lots of questions about corner cases. Here's the text from a section on rsIDs that I added into the appendix of my thesis; things have probably changed since then, but I expect the broad strokes are still correct:

Mike Feolo gave a talk at the [Working With The HapMap] Cambridge course on the dbSNP database. He talked about how new Single Nucleotide Polymorphisms (SNPs) were submitted, compared against currently existing SNPs, and merged if they were already in the database. The database RefSeq(rs) numbers are updated about twice a year, and whenever there is a new build of the human genome sequence. As a result of these builds, some rs numbers are merged to account for increased knowledge of mutations, a process that I was not aware of. This has an impact on my research, because I have received data from an Illumina 317k SNPchip, which I compare to the HapMap dataset. The rs numbers from the Illumina chip are for a specific dbSNP build, which may not be fully consistent with the most recent build (from which the HapMap data is retrieved). Mike has recommended that both the build number and the rs number are used when referring to a specific mutation. Sharma Buch pointed out that it was possible to ask Illumina for information about the merged rs numbers for their chips, which would allow datasets from different builds to be used together.

I think rsIDs are chr-pos specific, which is why multiple ALT alleles can exist in the same entry.

Thanks for all of this from all 3 of you. It's gradually helping me to form a picture. If anyone feels like reading and commenting on https://annovar.openbioinformatics.org/en/latest/articles/dbSNP/ it would be really useful. Is this guy an authority on this or not? His point that rsid means different thing to different people is revealing and if correct could prevent much head scratching and wasted hours reading and pondering this. It implies that one just has to get used to how rsid is used without trying to understand it too much (a sort of quantum mechanics like detachment...)

@Gringer your update is very interesting thanks.

Generally, consumer WGS companies do not deliver annotated VCF's. Annotations include rsIDs, gene names or similar area designations, and others.  But this can be done using the NIH dbSNP VCF files and the bcftools annotate command.  The annotation files at UCSC are modified for their GenomeBrowser and are often filtered and cleaned up.  While better and more consistent, they are usually not the latest release available from NIH.
Just make sure to match the build model and sequence (chromosome) naming convention between your annotation file and the original VCF. Dante usually delivers using the hs37d5 reference model. So use the .25 dbSNP VCFs files.  But dbSNP VCFs use accession names whereas Dante's hs37d5 reference model uses GRCh numeric ones.  So you will have to rename sequence names in one file or the other to make it work. Also, hopefully, the vendor has left aligned the InDels as that is what the dbSNP file has them defined as and what bcftools expects as well. Dante traditionally uses the DRAGEN pipeline with their Illumina sequencers; which does this. But they have recently switched to MGI DNB sequencers and the megaBOLT pipeline for customers outside the USA. We are not yet sure how they are processing InDels. Finally, Dante has historically delivered a separate SNP and InDel file.  You will often want to merge them or annotate them both separately then merge. bcftools can do this.
The dbSNP VCF file has well over 1 billion rsIDs now defined. And many fields of additional information that can be annotated in.  Often, if you just want to add a simple subset of the data, you can process the dbSNP VCF file to a simpler tab-separated (TSV, tab) file of entries. It can cut down the annotation file from 25 GB to 5 GB or less. See the bcftools annotation page for hints on using that tool to do this as well.
Personally, I find it easier to use online tools or services that have been tweaked.  Especially if they have nice browsers.  Sequencing dot com, usegalaxy.org and iobio.io are a few that come to mind.  Not all allow the export of the annotated VCF.
This and similar topics of utilizing your results are heavily discussed in the consumerWGS Facebook group; if you also wish to participate there.

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


text:
STACKS ref_map.pl no BAM records
I have 133 aligned and sorted BAM files that I am trying to run through the ref_map.pl script from the STACKS pipeline. My script looks as follows:
module load stacks

ref_map.pl -T 15 --popmap B2_filenames.txt --samples /home/heather/scratch/03_sort_and_map/filtered_and_sorted/B2/ --out-path /home/heather/scratch/ref_map_results/full_data

echo "ref_map.pl for B2 dataset completed"

I keep running into the issue where gstacks is unable to read the BAM files. This is what the log file looks like:
  Configuration for this run:
  Input mode: reference-based
  Population map: 'B2_filenames.txt'
  Input files: 133, e.g. '/home/heather/scratch/03_sort_and_map/filtered_and_sorted/B2/1401-47750.sort_filt.bam'
  Output to: '/home/heather/scratch/ref_map_results/full_data/'
  Model: marukilow (var_alpha: 0.01, gt_alpha: 0.05)

Reading BAM headers...
Processing all loci...
Error: No BAM records.
Error: (At the 0th record in file '/home/heather/scratch/03_sort_and_map/filtered_and_sorted/B2/1401-47754_B2.sort_filt.bam'.)
Aborted.

ref_map.pl: Aborted because the last command failed (1).
Last command executed by ref_map.pl was:
  /cvmfs/soft.computecanada.ca/easybuild/software/2020/avx2/Core/stacks/2.64/bin/gstacks -I /home/heather/scratch/03_sort_and_map/filtered_and_sorted/B2/ -M B2_filenames.txt -O /home/heather/scratch/ref_map_results/full_data -t 15

I checked the BAM file with samtools - the headers look ok, and I don't get any output when I run samtools quickcheck, so I can't quite figure out what is causing the issue. I am following the same pipeline as I did with a previous batch of data and can't understand whats going wrong.

"the headers look ok, and I don't get any output when I run samtools quickcheck" - this would be true even if the files contained only headers, no? Do you see alignment info when you look beyond the headers?

I keep getting a similar error with the sorted and filtered bam file generated from samtools "Error: 'sample.bam': this is a SAM file (and BAM was specified).


text:
PROKKA setupdb / cmpress error
I am trying to run --setupdb in Prokka, but I'm getting this error that I can't seem to find the answer for.
When I run
prokka --setupdb
I get
[18:57:31] Determined cmpress version is 001001 from '# INFERNAL 1.1.4 (Dec 2020)'
[18:57:31] Pressing CM database: /home/carol/.local/lib/prokka/db/cm/Viruses
[18:57:31] Running: cmpress \/home\/carol\/\.local\/lib\/prokka\/db\/cm\/Viruses

Error: File existence/permissions problem in trying to open CM file /home/carol/.local/lib/prokka/db/cm/Viruses.
CM file /home/carol/.local/lib/prokka/db/cm/Viruses not found (nor an .i1m binary of it)

[18:57:31] Could not run command: cmpress \/home\/carol\/\.local\/lib\/prokka\/db\/cm\/Viruses

Does anyone know how I can fix it?

This is a known problem with Prokka and was previously answered here https://bioinformatics.stackexchange.com/questions/20620/barrnap-bacterial-rrna-predictor-script-permission-denied-error-when-running-pro/20621#20621
Please let us know if this is okay, then we can close the question.


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

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

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

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

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

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

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

    print "##gff-version 3"
    next
}

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

from pyfaidx import Fasta

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

    return parser

@dataclass
class Feature:

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

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

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

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

def main():

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

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

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

        fasta_files = dict()
        seq_ids = dict()

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

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

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

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

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

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

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

            out.write(str(feature))

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

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

if __name__ == '__main__':
    main()

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

positional arguments:
  FILE                  The ABRicate TSV file

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

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

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

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


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

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

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

This sort of operation. Does that make sense?


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

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

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

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

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

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

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

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

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

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


text:
Manta Evidence Folder Is Empty
I am looking at Manta on Github and I ran a bit of code:
The following variables were defined in the code:

input_dir is the file path to a folder with less than 10 .bams as per the limitations of Manta
output_dir is the path to the output
MANTA_INSTALL_PATH is the path to the Manta installation
reference is the path to the reference fasta

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

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

# Now run each `runWorkflow.py`

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

The run worked and each sample had a folder with the following files:
results  runWorkflow.py  runWorkflow.py.config.pickle  workflow.error.log.txt  workflow.exitcode.txt  workflow.warning.log.txt  workspace

Each workflow.error.log.txt had nothing in it implying the the run had been successful. Each results folder had 3 subfolders
evidence  stats  variants

The stats and variants folder each had outputs but the evidence folder had nothing. Is this normal? Should it have something in this folder? Why is this folder there in Manta if nothing is in it?

The documentation on the GitHub page you linked to has this to say about the evidence folder:

Generating evidence bams
Using the --generateEvidenceBam option, Manta can be configured to
generate bam files of evidence reads for SVs listed in the candidate
vcf file.
It is recommended to use this option together with the --region
option, so that the analysis is limited to relatively small genomic
regions for debugging purposes.
The evidence bam files are provided in
${MANTA_ANALYSIS_PATH}/results/evidence, with a naming format
evidence_*.*.bam. There is one such file for each input bam of the
analysis, containing evidence reads of the candidate SVs identified
from that input bam. Each read in an evidence bam keeps all
information from the original bam, and it contains also a customized
tag in the format:
ZM:Z:${MANTA_SV_ID_1}|${EVIDENCE_TYPE},${MANTA_SV_ID_2}|${EVIDENCE_TYPE}.
For example, ZM:Z:MantaINV:5:0:1:0:0:0|PR|SRM,MantaDEL:5:1:2:0:0:0|SR
One read can have more than one of the three evidence types: PR for
paired reads, SR for split reads, and SRM for split read mates. One
read can be evidence for multiple SVs, which are separated by commas
in the tag. Notice that the number of evidence reads for a particular
SV in the evidence bam files could be more than the evidence counts
(PR and SR) in the final vcf files. This is because more stringent
criteria are applied for generating evidence counts in the final vcf
files.

So it looks to me like you need the --generateEvidenceBam option to generate evidence, but that the directory evidence is created regardless of whether you use the option. Try running the analysis again, but using --generateEvidenceBam this time.

Thank you I checked the documentation again this seems to be the case. It is just an empty folder until --generateEvidenceBam is specified.


text:
Singularity and Illumina's Nirvana
I am trying to using Illumina's Nirvana however I can't seem to get it to work despite having the right Singularity path.

I pulled Nirvana down using

singularity pull docker://annotation/nirvana:3.14

I found the path to the .dll using the following command:

singularity exec nirvana_3.14.sif /bin/bash -c "cd /opt/nirvana/ && ls"

Tried to run the downloader using:

## Define path to singularity image:
sif_container="/path/to/nirvana_3.14.sif"

## Path in the singularity file:
sif_path="/opt/nirvana/Downloader.dll"

## Data to be placed in
Data="path/to/Nirvana/GRCh38/"

#_____________________________________________________________________________________________________________________________________

singularity exec "$sif_container" dotnet "$sif_path" --ga GRCh38 --out "$Data"

yet it keeps saying ERROR: The top-level output directory directory was not specified. Please use the --out parameter.. I have checked the paths and they are correct and I have specified the data path yet it does not seem to recognize it? I have even moved around the container and the path variables to no avail?

From your previous questions, it is apparent that your HPC has an odd singularity setup. You could ask you sysadmin about how to navigate it.
EDIT
After working with OP on chat (and based on previous experience with them), it looks like their sysadmin has mount hostfs enabled in the singularity.conf file, which messes with every container every time. I guess it's easier than letting people bind their directories of choice but anyway, the way to circumvent that is to use --no-mount hostfs
@Indira, always run singularity [exec|shell|run] --no-mount hostfs instead of plain singularity [exec|shell|run]. All containers will then run well and if they error out, it's probably because either the container internally encountered an error or it was unable to access an external resource. What you're running into now is a mix that obscures the real problems and creates a no-win scenario for you.


text:
Singularity, AnnotSV and HPO
I am trying to use AnnotSV to annotate some Manta outputs in the following code with the HPO term changed:
Where:

input_dir is the input directory
singularity_container is the path to the AnnotSV singularity container

# Loop through the files in the input directory
for file in "$input_dir"/*.vcf; do
    # Get the base name of the file without extension
    filename=$(basename "$file" .vcf)

    # Create the output subdirectory within the input directory
    output_dir="$input_dir/AnnotSV_Batch_Run"
mkdir -p "$output_dir"

    # Set the output file path with "_AnnotSV" appended
    output_file="$output_dir/${filename}_AnnotSV_"

    # Run AnnotSV via Singularity with the base name as output prefix
    singularity exec "$singularity_container" AnnotSV -SVinputFile "$file" -genomeBuild GRCh38 --hpo "HP:0000001"  -outputFile "$output_file"
done

The HPO term I am using is in the jax.org site but when I run it - I get the most peculiar error  option not known. but the HPO term is specified in the AnnotSV manual on page 49 (according to the number on the bottom). As for trying I tried with --HPOterm and hpoterm I tried with and without the : and with and without quotes as well. I am not sure what the issue is but I cannot seem to find the solution. Has anyone perchance used this software before with HPO terms as it seems not to be accepting the input I give it?

The manual says -hpo and you're using --hpo. Could that be the problem?

Now I am getting "/usr/local/share/AnnotSV/Annotations_Human/Genes/GRCh38/refGene.txt.gz" file doesn't exist
Please check your install - Exit with error.?

Try adding a --no-mount flag immediately after singularity exec in the command.

The --no-mount gives the error FATAL:   "--no-mount": executable file not found in $PATH - I also tried with usr/local/bin/AnnotSV and it did not work

Oops, my bad. I did not read the documentation properly, sorry. Try singularity exec --no-mount hostfs AnnotSV

This error is persistent and I am not sure how to get around it.

Is this the page from which you're getting the singularity image link: https://lbgi.fr/AnnotSV/downloads

Try downloading the annotations file, uncompressing it and providing it as an input to the -annotationsDir parameter. The base image, which is ~400 MB seems to contain only the executable and all the images I see on docker hub are at least 4 gigs, which makes me think they include some sort of resource file that the base image does not have.

Yes the singularity is from the link

I believe I am doing this incorrectly as I give it the folder of the phenotype_to_genes.txt that I got from HPO but it seems like it still gives command not found for a lot this is only part of the error ```/usr/local/etc/AnnotSV/configfile: line 9: -annotationMode:: command not found
/usr/local/etc/AnnotSV/configfile: line 11: -benignAF:: command not found
/usr/local/etc/AnnotSV/configfile: line 13: -candidateGenesFiltering:: command not found
/usr/local/etc/AnnotSV/configfile: line 15: -genomeBuild:: command not found

The folder I gave was under the annotationsDir flag but the error persists.

Something is interfering with the shell parsing your command. Add an echo in front of the singularity command, then add a set -x before the for loop begins. Run the script and paste the output on pastebin, then link to the pastebin content here.

What is pastebin?

Did you Google it?

I did - it seems that AnnotSV doesn't work well with Singularity so this could not be done.

It seems that currently AnnotSV does not support HPO terms usage an Singularity.


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.


text:
GTF To Use With Broad Institute Fasta In STAR
To make a STAR reference genome used for alignment you need a genome fasta (the resources fasta provided by the Broad Institute in my case) and a GTF fasta. I tried to use the Ensembl human genome fasta fro hg38 but that does not work. I have looked in the GATK resources package but there is no .gtf file there. I really don't know what to do in this case? Is there anyone who can provide some guidance as to what .gtf file I should use?

If you are looking for a genome and its corresponding annotations, you can use the following, it is from nf-core (link to the original doc):
#!/bin/bash

VERSION=108
wget -L ftp://ftp.ensembl.org/pub/release-$VERSION/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa.gz
wget -L ftp://ftp.ensembl.org/pub/release-$VERSION/gtf/homo_sapiens/Homo_sapiens.GRCh38.$VERSION.gtf.gz

BTW the most up-to-date version as of today is 109.

When I use the Ensembl it doesn't work saying Fatal INPUT FILE error, no valid exon lines in the GTF file: and then the path to the Ensembl gtf.

That makes no sense, Ensembl's GTF is the most comprehensive GTF there is.

Thank you @gringer

I've had success in the past generating STAR indexes using the Gencode GTF files available from here:
https://www.gencodegenes.org/human/

Comprehensive gene annotation (ALL regions)
Genome sequence (ALL)

If you make sure that you're getting both the genome and the GTF files from the same source, it's more likely that they will work well together.
Note that the ALL annotation includes assembly patches and alternate haplotypes. The alternate haplotypes mean that population variation is captured, which usually increases the alignment rate of reads.
Unfortunately, most bioinformatics workflows don't properly handle alternate assemblies. As one example, reads aligning to those areas may end up classified as mapping to multiple genomic locations, even though they actually aren't. So it's a tradeoff of a more accurate alignment vs more hassle dealing with those alignments downstream.


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 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:
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 is the galaxy results file path?
I want to wrap galaxy into a docker image, and I want to delete files that are older than a certain date on a regular basis, so I want to know how to generate the files and where to store them.


text:
How to use Gblocks for trimming single-copy gene sequences?
I have run OrthoFinder on a set of 11 genomes and got the results in a folder. From the output folder, I saved all the single-copy gene sequences in a single FASTA file and now wish to remove all gaps using Gblocks before using them for constructing a phylogenomic tree. But I don't know the command line used for running Gblocks. Any help will be appreciated!

I assume you've aligned these sequences? I don't know GBlocks but I know trimal and will do the same thing (or do it better). Could you describe what you mean by 'removing gaps'? 1. Are you looking to trim the 5' and 3' ends? 2. Wanting to remove 'gap [indel] only columns'? 3. May be wanting to remove the any position which contains an indel, i.e. the entire site is removed? The last option is the best approach.

Did you look for the gblocks docs? You know it's >20 years old that tool!

Hey @M__ ! I am still learning the procedures for doing phylogenetic analysis but I think I need to remove ambiguous sequences that hamper the accuracy of a phylogeny. You can check this paper which has cited Gblocks.                                        [link] (https://doi.org/10.1016/j.syapm.2020.126166)          If you feel trimal is better could you please explain how to use it? Thanks!

@Chris_Rands, I did come across the Gblocks documentation but it did not work for me. I think it has not been updated in a while. Since my supervisor shared the paper which used Gblocks, I wanted to learn it. I agree that the latest tool would be better.

Basically this is 16S data. If its proteins there's a different strategy - even if the data set needed is nucleotides (separate question). I would simply use trimal which is available on conda
trimal -in <inputfile> -out <outputfile> -nogaps

The input file is aligned fasta here and will output as aligned fasta.
Any site (alignment position) with an indel - will be removed. This is definitely preferred if performing an iqtree analysis, because of the way it handles indels. This will also clean up poorly aligned, alignments.
e.g.
---ATG
---A-G
--TATG

This will result in...
AG
AG
AG

Ultimately for alignments - there's no magic bullet - and manual editing and inspection is expected e.g. using seqotron. There are very good reasons for this concerning circular mathematics. However, modern aligners such as muscle5 are extremely good.
Trimal is available via conda
conda install -c bioconda trimal

In terms of using GBlock I'd need to know whether it was simply deleting indel sites, deleting indel sites at >50% frequency or whether it was performing a genetic diversity calculation per site. I don't think it is doing the latter. It is better to perform diversity saturation calculation explicitly rather than using an unknown automate because thats kinda important. If its deleting sites with excessive >50% indels - honestly its better to do a blockwise deletion rather than do something half-half. Within SSU there is a reason might be being done because there can be excessive indels at deep taxonomic levels with this data type. Again using iqtree just do it 'blockwise' which is the trimal output.
Note I've not read the original GBlock paper, but if you summarise its application I could consider it against trimal. However, I definitely prefer trimal considering the downstream application.

Thank you again, @M__ . Will definitely go ahead with what you have suggested. Although, I must admit, being an absolute amateur in bioinformatics I don't follow some of the things you have so sincerely explained. Forgive me if I am being imprudent, do you have any way that I could learn all this under your guidance? I wish to be thorough.. thanks again!

If you need the theory you can always ask a separate question. The theory is very detailed, why certain things and done and otherwise not.


text:
Error with eggNOG-mapper jobs
I recently tried using eggNOG-mapper for pan-genome and core genome analysis of files from OrthoFinder, specifically the 'Orthogroup_Sequences' which contains a huge list of protein files in fasta format. As instructed, the folder was compressed into a tar.gz zip file and uploaded. If the upload was successful, an email is sent to your mail ID where you can check the details of your job and decide whether to START the job or not. Every time I START the job, it fails within a few minutes. Initially I presumed that it was due to the file size and hence reduced the number of sequences submitted, but in vain. Can anyone help me?
eggNOG-mapper : http://eggnog-mapper.embl.de/

Which version of eggNOG-mapper is being used here?

Hi @M__ I think it's v2. That's what the citation says at the bottom of the page. it's probably v2.1.12. Do you think the command line version would be better?

There is a known bug between eggNOG-mapper2 and eggNOG6. This is described here:
https://github.com/eggnogdb/eggnog-mapper/issues/428
In your case, I'm not sure on what strategy you are deploying after OrthoFinder, but there's a chance it's not solvable at present and the best place for information is the eggnog GitHub.

Thank you @M__ . I have posted my issue on the github page for failing eggNOG-mapper jobs. I have also tried to perform the same function using COGclassifer after concatenating all the .fa sequences in the Orthogroup_sequences file. Hopefully it yields reliable results. If you have any other suggestions please let me know.

@K_081 Thanks it's a positive step forward. eggNOG is very, very good. I'm not be looking for a solution outside this. I'll keep an eye out for your GitHub post and I will be looking around for a way forward and post back here once done (timeline several weeks).

Hi @M__, I just got a reply on the github page from one of the developers I guess, and he has found an unwanted line in my file.  It is the file created by OrthoFinder itself so I don't know how it appeared there. I'll have to check all the files before uploading them to eggNOG-mapper.

I'd like to answer my own question as I have found the solution for this and thought it might help anyone else trying to use eggNOG-mapper.
The prerequisite for uploading the input file with the required collection of .faa files, is that, the file extension has to end with .gz. The site mentions that too. I had made the mistake of uploading a .tar.gz which is what my system produced as a compressed file of my folder of sequences. Although the upload was successful, the runs failed as the file had some unwanted lines owing to the .tar.gz format.
I randomly uploaded a concatenated .faa file which worked well and got the results too! The trick is that, the input file has to be strictly a .gz and NOT .tar.gz or it could also be a .faa file with all the sequences.

Thanks for this that is very helpful.


text:
How to access specific SRA sequences within an SRP file?
How do you get a specific 100 SRA sequences from a Bioproject which contains 400 runs, when you don't need all of them?

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

Welcome @kirtikulanthaivel, the answer is very likely to involve the SRA toolkit (freeware), if it is doable. How to apply this is more challenging because the question isn't clear. Could you supply a reproducible example? It does not need to be the data set(s) you are interested in.

I actually wanted to know how a specific number of SRR runs, for example 100 of them from a total of 400 in an accession list can be downloaded at a time. Problem solved, thanks to a friend. Thank you all anyway.


text:
Is it normal to have low mapping to bacterial genome in a total RNA sample of plant root-colonized bacterial cells?
I had isolated the total RNA from a sample of plant roots with colonized bacteria by excising the roots (harboring the bacterial cells) and crushing them into a powder using liquid nitrogen. The RNA isolation procedure was done following the kit's instructions specific for bacterial RNA.
The RNA was sequenced and the raw reads obtained were checked for mapping percentage with the bacterial genome's GTF file using bowtie2, as I need to study the expression of the bacterial genes. While the mapping was close to 80% for samples coming from the planktonic state, it was quite less - 0.4 to 0.56%, for all the samples coming from root colonized cells, obviously because a major part of it is comprised of plant RNA. Firstly, is this acceptable? Can I use these reads for further analysis?
Secondly, I use DESeq2 to analyze the data for upregulated and downregulated genes, comparing the colonized cells to the planktonic ones. Since the package includes a normalization step, is it safe to assume that although the mapping is less it does not influence the analysis significantly?


text:
Does using a microarray chip that matches the chip used in training lead to higher explained variance of a polygenic score?
I'm trying to replicate an existing polygenic score (i.e. test the accuracy in a new sample), and want to know if matching the original study's microarray chip will improve the accuracy (that is, variance explained in a new test set) of the polygenic score.
The labs I have access to have different multiple microarray chips to choose from (e.g., GSA, GDA, UK Biobank Axiom Array, UK BiLEVE Axiom array, GeneChip 2.0 array, Infinium Core-24 Kit, etc.). Will using a chip that is different than the chip used in the study that produced the polygenic score impact the explained variance (accuracy) of the polygenic score? If so, by how much?
Example: a study on heart disease runs a GWAS on 900,000 individuals who are genotyped with the UK BiLEVE Axiom array. This study uses these summary statistics to produce weights for a polygenic score. Lab A, B, and C are conducting polygenic score validation/replication studies, and they each have a random sample of 10,000 individuals drawn from the general population, and they all use the same summary statistics from the original study. After imputing the genotype data of the 10,000 people, they calculate the polygenic score for each person, and then evaluate how much variance in the phenotype they are able to explain. The only difference is that Lab A genotypes their 10,000 people with the UK BiLEVE Axiom array, whereas lab B uses the Illumina GSA, and lab C uses the Illumina GDA.
Do we have any reason to suspect that there will be a substantial difference in variance explained (r^2) by the polygenic score between labs? Are there any studies which have analyzed something that could be informative to this question?
I've been unable to find research on this question. This should not be dependent on which polygenic score is being studied; there is no reason to specify a particular phenotype.
Note: it is standard to impute both in the original study and the replication.

Genotyping can be validated by whole-genome sequencing, and this has been done on the International Genome Sample Resource. You should be less concerned about the reproducibility of genotyping results, and more concerned about the correctness or validity of genotyping results. Companies will often revisit old SNPchips to mark invalid spots, update rs numbers, and add new genotyping probes that are better on a population-wide scale. In that sense, choosing a SNPchip platform that has an established track record of QC and information about version updates would be a better idea than trying to match a previous study as precisely as possible.
However, with regards to replicating existing studies, a bigger concern for you should be the variation in sampled populations, which plays a big role in the outcome of a polygenic test.
This variation is so extreme that directions of associations can flip for the same variant, especially for low-frequency variants with extremely low p-values (i.e. the ones that often show up as manhattan plot peaks).
The issues around population and study replication are explained well in the following paper from 2001:
https://doi.org/10.1006/tpbi.2001.1543
This is not something that can be fixed by attempting to replicate studies in more populations. As Cecile Janssens says, larger sample size only affects random error, not systematic error.
https://twitter.com/cecilejanssens/status/1119235476567801856


text:
How do I infer sex from a VCF file?
I have human samples VCF files that I want to infer the sex of. I do not have sequencing information or information about coverage. What's the best way to go about this? I'm looking for a robust and reliable method. Tools are okay, and conceptual or algorithmic ideas about what to do are also okay.

try somalier: https://github.com/brentp/somalier
https://genomemedicine.biomedcentral.com/articles/10.1186/s13073-020-00761-2 :

The sites files include informative variants on the X and Y chromosomes so that Somalier can also estimate a sample’s sex from the genotypes.


text:
StringTie discovers less novel isoforms when reference annotation is provided
I am analyzing PacBio IsoSeq data and I am using StringTie to assemble transcripts. I noticed the difference in the output files when I add a reference genome annotation (-G option in StringTie). In general I have less transcripts assembled when I give as imput reference genome annotation and I was wondering why. I know that the annotation of the genome is not the best (Borelian orangutang) and I would rather trust the data rather than the annotation provided.
I am attaching an example for one of the genes (IGV): the first track is reference annotation; the second is StringTie output when no reference annotation is provided; the third is StringTie output when reference annotation is provided (in the third case two novel assembled isoforms are missing)

Not an answer. For novel transcripts, you may consider IsoQuant. Its Nature Biotech paper argued that it does particularly better than StringTie on novel transcripts.

Not very familiar with what StringTie is doing, but in general reference annotations are helpful in providing a prior for what the tool expects.
I would suggest looking at the manual and using some of the provided options and outputs to probe the coverage of the various transcript species (I believe these are in the "attributes" section of the GTF output).
If the transcripts that are missing when you use -G have a coverage of just a couple of reads, they are likely to be less interesting. On the other hand, if they have relatively high coverage, then you have more reason to be skeptical of the existing annotation.
My guess is that the tool is more heavily weighting existing transcripts based on the annotation, and therefore excluding the extra transcripts as being below some threshold.
I'd recommend consulting the "recommended workflow" in the manual that supposedly allows for novel isoforms. Ensure that you are using this workflow and not the more streamlined workflow, which does not allow for novel isoforms.
See also the comment suggesting comparing to other IsoSeq-oriented tools.


text:
maxquant error: Unhandled Exception: System.IndexOutOfRangeException: Index was outside the bounds of the array
An error occurred while I was running the program maxquant (2.0.3):
Reading search engine results (SP)
Unhandled Exception:
System.Exception: Exception during execution of external process: 1228404
Unhandled Exception:
System.IndexOutOfRangeException: Index was outside the bounds of the array.
at MsLib.Search.AndromedaPeptide.FromString (System.String line, BaseLibS.Ms.Data.Protein.ProteinSet proteinSet) [0x000a9] in :0
at MsLib.Search.AndromedaQuerySp.Add (System.Int32 count, System.String line, BaseLibS.Ms.Data.Protein.ProteinSet proteinSet, System.Int32 minLength) [0x00001] in :0
[ERROR] FATAL UNHANDLED EXCEPTION: System.IndexOutOfRangeException: Index was outside the bounds of the array.
at MsLib.Search.AndromedaPeptide.FromString (System.String line, BaseLibS.Ms.Data.Protein.ProteinSet proteinSet) [0x000a9] in :0
how can i fix it, thank your very much!

Thank you for your reply. I only changed the modification and added the fasta  raw file  in the configuration file, Other parameters were default.

Dear @jingzhao thanks for posting back. Could you kindly "accept" your answer (green tick button) - it does help with our site stats

****this was the xml file that i configured to maxquant. i changed modification and added the fastq and raw file. Other parameters were default. however, 'Unhandled Exception: System.IndexOutOfRangeException: Index was outside the bounds of the array' was reported when I run the command 'maquant mq.xml'. but, When I only kept the N modification, the program ran successfully, but a small amount of peptides were identified.


text:
Generating 3D cartesian coordinates of a peptide and computing van der Waal interactions (1-4 atoms) for a given set of dihedral angles?
I wanted to generate 3D Cartesian coordinates of a peptide sequence for fixed bond lengths and bond angles for a given set of backbone dihedral angles and side dihedral angles. I also wanted to compute non-bonded van der Waal interactions (1-4 atoms) for each of the peptide generated from the dihedral angles. I initially tried frag builder module in python although it works fine. There are few issues,

I could not fix bond length and bond angle of my own choice.

It generates PDB file and save it inside the hard disk. So I need to read it again and parse it to compute non-bonded van der Waal interactions (1-4 atom). It consumes lot of time

Also I find fragbuilder is not that fast in generating peptides.

Reading , parsing and computing the interactions takes around 30 mins to complete the task for 200,000 3D ensemble of structures of given peptide sequence of length 18.

Kindly suggest me a better module in python which does the job faster and also give reference to it.

Don't ask multiple questions in a single post. Focus on solving one problem at a time.

@M__    Yes Fragbuilder has many shortcomings as pointed out in the questing itself. So I am looking for better alternative. Yes, I am limited to python. Is there any other modules which could handle this inside python ?  Yet, I want to know whether there is any other method outside python as well. Thanks in advance !

Thanks @vigneshwarankannan updating the question would be cool.

The PDB module in Biopython is worth checking out for obtaining angles, distances, torsion angles. It is described from p208 at http://biopython.org/DIST/docs/tutorial/Tutorial.pdf. It doesn't do van der Waal interactions.
For example,
from Bio.PDB.PDBParser import PDBParser
from Bio.PDB.Chain import Chain
from Bio.PDB.internal_coords import *
from Bio.PDB.PICIO import write_PIC, read_PIC, read_PIC_seq
from Bio.PDB.ic_rebuild import write_PDB, IC_duplicate, structure_rebuild_test
from Bio.PDB.SCADIO import write_SCAD
from Bio.PDB.PDBIO import PDBIO

parser = PDBParser(PERMISSIVE=1)
myProtein = parser.get_structure("1a8o", "1A8O.pdb")
myChain = myProtein[0]["A"]
IC_Chain.MaxPeptideBond = 4.0
myChain.internal_coord = None
myChain.atom_to_internal_coordinates(verbose=True)
r2 = myChain.child_list[1]
r2ic = r2.internal_coord
r2chi1 = r2ic.get_angle("chi1")

I wanted to compute the 3D structures of the peptides (i.e) x,y,z coordinates of the atoms for the given set of input Ramachandran angles for a fixed bond lengths and bond angles. I don't require to compute the Ramachandran angles of a given PDB file.

@vigneshwarankannan thanks for the feedback. What you want is simply to extend the Biopython module with appropriate code for the trigonometric calculation, i.e. angles and distances. It's not my area - but it would be doable.


text:
Is there any python module to assign secondary structures which includes both Polyprolene region as well as individual strand as well?
I wanted to assign secondary structure to set of pdb files. However, DSSP and STRIDE cannot assign for single peptides if it is falling inside single long beta strand or polyprolene region as it does not have any hydrogen bond.DSSP and STRIDE work fine as long as the peptides form hydrogen bonds.
Can someone suggest me a way to assign secondary structure for peptides which also fall in the region of PPII and beta strand region ?