text:
Do proper .fcs flow-cytometry files specify the gating hierarchy, so one may calculate cell-counts without specifying it themselves?
Do proper .fcs flow-cytometry files already have the gating hierarchy, so one may calculate cell-counts without extra work?
https://github.com/whitews/FlowKit/
# setup some file paths for our data
base_dir = "data/8_color_data_set"

sample_path = os.path.join(base_dir, "fcs_files")
wsp_path = os.path.join(base_dir, "8_color_ICS.wsp")
session = fk.Session(sample_path)
session.import_flowjo_workspace(wsp_path)
# that's it, let's see what sample groups were loaded
session.get_sample_groups()
['default', 'All Samples', 'DEN']
# I happen to know all the "real" analysis is within the "DEN" group
sample_group = 'DEN'
# and the sample that are included in the group
session.get_group_sample_ids(sample_group)
['101_DEN084Y5_15_E03_009_clean.fcs',
 '101_DEN084Y5_15_E01_008_clean.fcs',
 '101_DEN084Y5_15_E05_010_clean.fcs']
# so what is the gating hierarchy for this group
print(session.get_gate_hierarchy(sample_group))
root
╰── Time
    ╰── Singlets
        ╰── aAmine-
            ╰── CD3+
                ├── CD4+
                │   ├── CD107a+
                │   ├── IFNg+
                │   ├── IL2+
                │   ╰── TNFa+
                ╰── CD8+
                    ├── CD107a+
                    ├── IFNg+
                    ├── IL2+
                    ╰── TNFa+
# looks good, let's go ahead and analyze them (using verbose mode to see each gate as it's processed)
session.analyze_samples(sample_group, verbose=True)
#### Processing gates for 3 samples (multiprocessing is enabled - 3 cpus) ####
101_DEN084Y5_15_E03_009_clean.fcs: processing gate Time
101_DEN084Y5_15_E03_009_clean.fcs: processing gate Singlets
101_DEN084Y5_15_E03_009_clean.fcs: processing gate aAmine-
101_DEN084Y5_15_E01_008_clean.fcs: processing gate Time
101_DEN084Y5_15_E01_008_clean.fcs: processing gate Singlets
101_DEN084Y5_15_E01_008_clean.fcs: processing gate aAmine-
101_DEN084Y5_15_E05_010_clean.fcs: processing gate Time
101_DEN084Y5_15_E03_009_clean.fcs: processing gate CD3+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate Singlets
101_DEN084Y5_15_E05_010_clean.fcs: processing gate aAmine-
101_DEN084Y5_15_E03_009_clean.fcs: processing gate CD4+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate CD3+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate CD4+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate CD3+
101_DEN084Y5_15_E03_009_clean.fcs: processing gate CD8+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate CD4+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate CD8+
101_DEN084Y5_15_E03_009_clean.fcs: processing gate CD107a+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate CD8+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate CD107a+
101_DEN084Y5_15_E03_009_clean.fcs: processing gate IFNg+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate CD107a+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate IFNg+
101_DEN084Y5_15_E03_009_clean.fcs: processing gate IL2+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate IFNg+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate IL2+
101_DEN084Y5_15_E03_009_clean.fcs: processing gate TNFa+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate IL2+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate TNFa+
101_DEN084Y5_15_E03_009_clean.fcs: processing gate CD107a+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate TNFa+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate CD107a+
101_DEN084Y5_15_E03_009_clean.fcs: processing gate IFNg+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate CD107a+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate IFNg+
101_DEN084Y5_15_E03_009_clean.fcs: processing gate IL2+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate IFNg+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate IL2+
101_DEN084Y5_15_E03_009_clean.fcs: processing gate TNFa+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate IL2+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate TNFa+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate TNFa+
# and a look a the results
session.get_group_report(sample_group)
sample  gate_path   gate_name   gate_type   quadrant_parent parent  count   absolute_percent    relative_percent    level
0   101_DEN084Y5_15_E03_009_clean.fcs   (root,) Time    RectangleGate   None    None    283968  99.999648   99.999648   1
1   101_DEN084Y5_15_E03_009_clean.fcs   (root, Time)    Singlets    PolygonGate None    Time    236780  83.382341   83.382635   2
2   101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets)  aAmine- PolygonGate None    Singlets    161823  56.986150   68.343188   3
3   101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-) CD3+    PolygonGate None    aAmine- 132200  46.554377   81.694197   4
4   101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-, CD3+)   CD4+    PolygonGate None    CD3+    81855   28.825330   61.917549   5
5   101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-, CD3+)   CD8+    PolygonGate None    CD3+    46965   16.538777   35.525719   5
6   101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) CD107a+ RectangleGate   None    CD4+    59  0.020777    0.072079    6
10  101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) CD107a+ RectangleGate   None    CD8+    814 0.286651    1.733206    6
7   101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) IFNg+   RectangleGate   None    CD4+    598 0.210586    0.730560    6
11  101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) IFNg+   RectangleGate   None    CD8+    1224    0.431033    2.606196    6
8   101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) IL2+    RectangleGate   None    CD4+    64  0.022538    0.078187    6
12  101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) IL2+    RectangleGate   None    CD8+    127 0.044723    0.270414    6
9   101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) TNFa+   RectangleGate   None    CD4+    288 0.101420    0.351842    6
13  101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) TNFa+   RectangleGate   None    CD8+    821 0.289116    1.748110    6
0   101_DEN084Y5_15_E01_008_clean.fcs   (root,) Time    RectangleGate   None    None    290166  99.997932   99.997932   1
1   101_DEN084Y5_15_E01_008_clean.fcs   (root, Time)    Singlets    PolygonGate None    Time    239001  82.365287   82.366990   2
2   101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets)  aAmine- PolygonGate None    Singlets    164655  56.743931   68.893017   3
3   101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-) CD3+    PolygonGate None    aAmine- 133670  46.065782   81.181865   4
4   101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-, CD3+)   CD4+    PolygonGate None    CD3+    82484   28.425899   61.707189   5
5   101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-, CD3+)   CD8+    PolygonGate None    CD3+    47165   16.254153   35.284656   5
6   101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) CD107a+ RectangleGate   None    CD4+    68  0.023434    0.082440    6
10  101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) CD107a+ RectangleGate   None    CD8+    73  0.025157    0.154776    6
7   101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) IFNg+   RectangleGate   None    CD4+    4   0.001378    0.004849    6
11  101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) IFNg+   RectangleGate   None    CD8+    2   0.000689    0.004240    6
8   101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) IL2+    RectangleGate   None    CD4+    6   0.002068    0.007274    6
12  101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) IL2+    RectangleGate   None    CD8+    10  0.003446    0.021202    6
9   101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) TNFa+   RectangleGate   None    CD4+    21  0.007237    0.025459    6
13  101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) TNFa+   RectangleGate   None    CD8+    8   0.002757    0.016962    6
0   101_DEN084Y5_15_E05_010_clean.fcs   (root,) Time    RectangleGate   None    None    284846  99.844369   99.844369   1
1   101_DEN084Y5_15_E05_010_clean.fcs   (root, Time)    Singlets    PolygonGate None    Time    235331  82.488345   82.616923   2
2   101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets)  aAmine- PolygonGate None    Singlets    160955  56.418031   68.395154   3
3   101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-) CD3+    PolygonGate None    aAmine- 131305  46.025097   81.578702   4
4   101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-, CD3+)   CD4+    PolygonGate None    CD3+    81135   28.439483   61.791249   5
5   101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-, CD3+)   CD8+    PolygonGate None    CD3+    46546   16.315328   35.448764   5
6   101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) CD107a+ RectangleGate   None    CD4+    70  0.024536    0.086276    6
10  101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) CD107a+ RectangleGate   None    CD8+    306 0.107259    0.657414    6
7   101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) IFNg+   RectangleGate   None    CD4+    9   0.003155    0.011093    6
11  101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) IFNg+   RectangleGate   None    CD8+    348 0.121981    0.747647    6
8   101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) IL2+    RectangleGate   None    CD4+    8   0.002804    0.009860    6
12  101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) IL2+    RectangleGate   None    CD8+    65  0.022784    0.139647    6
9   101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) TNFa+   RectangleGate   None    CD4+    28  0.009815    0.034510    6
13  101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) TNFa+   RectangleGate   None    CD8+    336 0.117775    0.721867    6

The flow cytometry standard (FCS) file specification does not include gating information. There is a separate specification, GatingML, that does provide gating definitions, though few applications adhere strictly to that standard. As the author of FlowKit, we wanted to provide a tool for the Python bioinformatics community that reads and creates valid GatingML documents for FCS data.


text:
How does one convert "flowJO" flow-cytometry cell-count, etc, bead, etc data into .json or relational-database-management-system compatible files?
How does one convert "flowJO" flow-cytometry cell-count data into ".json" or relational-database-management-system compatible files?
The way flow-jo exports is in wide rather than long format. It also exports everything per subject/animal/lab-rat per specific timepoints rather than aggregating. I need the aggregated files to plot groups.
It's also not shaped in a "primary-key", column 1, column 2, column 3, or .json (NO_SQL) format.
Running statistical programs or databases is intractable with the export options Flowjo gives.
I am looking for some way to convert it, or maybe the problem lies in flowjo. It does not support the export formats I need.
I amk open to other flow-cytometry software that offers better integration with no-sql/sql 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.

If you are familiar with Python, I maintain the FlowKit library that can import FlowJo 10 workspaces and export the gate statistics as a Pandas DataFrame. From there it would be trivial to convert the data to JSON or CSV.
If you want to go that route, the Workspace class handles this functionality. There are tutorials to get you up to speed, the most relevant is this Jupyter notebook demonstrating loading a FlowJo .wsp file along with the related FCS files and extracting the gate analysis.

source: Stack Exchange

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

sort: votes


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

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

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


text:
Can I discard missing alleles?
I am converting a biallelic VCF into a SQL table and one of my tables will be something like:
SAMPLE | custom variant identifier | Genotype
-------|---------------------------|---------
0001   | 3456                      | ./.

My question is that for the diploid ./. or haploid . missing alleles, whether I can just drop them to save space. My thinking is that since the dot represents missing data, the ./. are just placeholders in the VCF but provide no information at all.
Is this correct? Should I retain the rows with ./.?
Showing the end of one line of the VCF (really long lines because of high N):
FORMAT          SampleX             SampleY
GT:AD:DP:GQ:PL  0/0:3,0:3:4:0,4,45  ./.:3,0:3:.:0,0,0

I thought to also show the result after bcftools norm -m to split to biallelic and queried per sample.
sample_id   chr@pos@ref@alt genotype
sample1 chr1@10158@ACCCT@A       0/0
sample2 chr1@10158@ACCCT@A       ./. <- does this provide any info or can I discard?
sample3 chr1@10158@ACCCT@A       0/0
sample4 chr1@10158@ACCCT@A       ./.
sample5 chr1@10158@ACCCT@A       0/0

Can you show us a few lines of the actual VCF file please? ./. means the site wasn't genotyped. It is common in gvcf files, for example but we'd have to see some data to judge correctly. Where are the "asterisks"? Do you mean .? Asterisks are *.

Thanks, I've updated the question. I guess what I'm asking is whether there is any information being encoded for the samples that weren't genotyped for that variant i.e. ./.

If you can, do not discard missing allele information.
Generally you often benefit from being able to distinguish between a position that you interrogated and genotyped as REF (reference) versus one that you interrogated but could not determine the genotype (e.g due to low coverage in the sequencing or probing issues in a microarray setup) and which usually is encoded as a missing (./. if ploidy = 2) genotype. A missing genotype tell us that the genotype at this position could be REF, but could also be some other base, we simply do not know.
Knowing this is critical and/or useful depending the downstream analysis. For example it is quite important and informative when merging datasets. Since formats like VCF or PLINK do not record non-variant positions, if you discarded missing genotypes prior to merging, how can you tell between genotypes that truly were REF in one dataset (fully invariant position) and some recorded missing genotypes at the same position but which you decided to discard? You will end up assuming REF status for everything not recorded and this will inevitably lead to batch effects in your combined dataset. For a better discussion see this blog post
On the other hand, if your relational database (your particular situation) for instance will not incorporate external databases in order to produce joint datasets, perhaps discarding missing genotypes is all fine.


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

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

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

Part 2

of each eigengen is around 30%.

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

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


text:
What are the columns in bedtools coverage output hist file?
I am using bedtools to caculate the coverage of my targets of my WES data, to later plot. But to plot, I need to know what to plot and what it is I am seeing. I have unsuccefully tried to find what exactly the columns are in the output of the bedtool coverage, and how to read them. Also, the documentation does not help me a lot, and I find it confusing. Does anyone know, or where to find this a bit more elaboratly explained?
I used the following code:
bedtools coverage -hist -abam input.bam -b input.bed > output >2 log
Correction note: the code above is actually wrong (thanks to @Alexlok) and should be:
bedtools coverage -hist -a input.bed -b input.bam > output 2> log
And the output is (first line):

Last 5 lines:

Thank you in advance

Make sure that you're not following a script written for bedtools version < 2.24: the inputs a and b are inverted in recent versions. So on a recent version of bedtools, if you want to count the read coverage of your features, you need: bedtools coverage -a input.bed -b input.bam. Your current command counts the number of features in b that intersect each read in a, and I'm not sure that what you want to do.

I am following workflows for exome coverage probally indeed using <2.24, since they were published a while ago... not finding any new versions. I want to see how much of the exome targets (bed file) are actually captured by my reads indeed. I use -abam (my aligned reads to the reference genome) and -b (my target bed file) option? So I am not actually using the -a option. Does it also count then? Thanks

If I understand correctly (I have only used bedtools a bit, no expert on that), on older versions (before 2.19.0) -a would have expected a bed file, so -abam made explicit that it was a bam. But on more recent versions, -a and -b automatically detect the input type, so -abam is synonym with -a (kept for compatibility). And since version 2.24 the options a and b are inverted, so you would need a -bbam (which does not exist since -b covers that use case). See the example in the page you linked.

Thank you so much for that. You are absolutly right and I got totally confused.

Found it, it is an extended .bed file format, so not specific to bedtools output: source bed file output.
With 4 added columns to the end from bedtools -hist option:
bedtools -hist option

Thanks for posting the answer, appreciated.


text:
Multiple values for memory resources shown when Snakemake is run
I am running a Snakemake (V7.14.0) pipeline on the cluster and provided resources in a resource config file like:
threads:
    tool1: 4
mem_gb:
    tool1: 12

Point to it in my rules like:
resources:
        mem_gb=config['mem_gb']['tool1']
threads:
        config['threads']['tool1']

And the snakemake cluster command like:
snakemake --configfiles \
    config/config_resources.yml \
    config/config_parameters.yml \
    --latency-wait 60 \
    --use-conda \
    --use-singularity \
    --printshellcmds \
    --jobname MyPipe{name}_{jobid} \
    --jobs 300 \
    --cluster \
    "bsub -W 200 -n {threads} -M {resources.mem_gb}G \
    -R 'rusage[mem={resources.mem_gb}G]' \
    -o run_logs/{name}_%J.out \
    -e run_logs/{name}_%J.err"

What I do not understand when I run the pipeline is that it shows this in the terminal:

threads: 4
resources: mem_mb=19442, disk_mb=19442, tmpdir=/tmp, mem_gb=12

So next to the one I provided (mem_gb), there are: mem_mb, and disk_mb. They are clearly not the same value. In the cluster job it looks as though the MEMLIMIT is indeed the one I give, but the MEM usage is different (higher). Is it because Snakemake does not understand the mem_gb and I specifically need to use that word?
Thanks in advance

You have defined mem_gb but not mem_mb in your resource configuration file. You can set them to same value by adding something like mem_mb=config['mem_gb']['tool1']. Additional information in their documentation.

Thanks for the response. So basically, Snakemake does not understand mem_gb. So to do this correctly, I need to actually *1000 the values. Because otherwise I get too low values (e.g. 6 mb instead of 6 gb).


text:
What should GC coverage bias plot of exome data look like?
I want to see if there is any GC coverage bias in some paired-end Illumina human exome data. Update This is for a variant calling project. I used the CollectMultipleMetrics from GATK4 toolkit with -PROGRAM CollectGcBiasMetrics. An example of the command I used is:
gatk CollectMultipleMetrics -I input.bam -O prefix_dir_name -R input.ref_genome -PROGRAM CollectGcBiasMetrics -PROGRAM CollectAlignmentSummaryMetrics -PROGRAM CollectInsertSizeMetrics -PROGRAM QualityScoreDistribution -PROGRAM MeanQualityByCycle -PROGRAM CollectBaseDistributionByCycle

Update This is my result for the QC coverage bias with one sample highlighted:

There seems to be a large bias. I would expect the peak to be around 50% for human. And those weird peaks around 90% are also puzzling to me. On the other hand, I am working with exome data so it might be logical to have a bias towards the higher values of GC% because you have reads on GC-rich exons?
Update I have another plot with the bins produced by the above mentioned tool (.gc_bias.pdf)

Update first part of useful file (gc_bias.detail_metrics) which shows the reads per bin, as suggested by @Maximilian Press. This I will plot.
## htsjdk.samtools.metrics.StringHeader
# CollectMultipleMetrics --INPUT WESP_run/results/aligned/RM8398_MF_S8_MD.sorted.bam --OUTPUT WESP_run/results/aligned/post_alignment_QC/multiple_metrics/RM8398_MF_S8 --PROGRAM CollectAlignmentSummaryMetrics --P
ROGRAM CollectBaseDistributionByCycle --PROGRAM CollectInsertSizeMetrics --PROGRAM MeanQualityByCycle --PROGRAM QualityScoreDistribution --PROGRAM CollectGcBiasMetrics --REFERENCE_SEQUENCE /mnt/scratch_dir/deute
koe/projects/HP/WESPipe/workflow/resources/reference_genomes/hg38/alt_plus/GCA_000001405.15_GRCh38_full_plus_hs38d1_analysis_set.fna --ASSUME_SORTED true --STOP_AFTER 0 --METRIC_ACCUMULATION_LEVEL ALL_READS --IN
CLUDE_UNPAIRED false --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 500000 --CREATE_INDEX false --CREATE_MD5_FILE false --GA4GH_CLIENT_SECRETS client_sec
rets.json --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false
## htsjdk.samtools.metrics.StringHeader
# Started on: Sat Dec 03 14:25:15 CET 2022

## METRICS CLASS        picard.analysis.GcBiasDetailMetrics
ACCUMULATION_LEVEL      READS_USED      GC      WINDOWS READ_STARTS     MEAN_BASE_QUALITY       NORMALIZED_COVERAGE     ERROR_BAR_WIDTH SAMPLE  LIBRARY READ_GROUP
All Reads       ALL     0       139942  1326    19      0.413544        0.011357
All Reads       ALL     1       102412  61      17      0.025996        0.003328
All Reads       ALL     2       118957  74      16      0.02715 0.003156
All Reads       ALL     3       152181  929     20      0.266429        0.008741
All Reads       ALL     4       167732  248     16      0.06453 0.004098
All Reads       ALL     5       173319  232     16      0.058421        0.003836
All Reads       ALL     6       196365  304     15      0.067567        0.003875
All Reads       ALL     7       211217  343     17      0.070875        0.003827
All Reads       ALL     8       241125  441     18      0.079822        0.003801
All Reads       ALL     9       270272  1068    18      0.172463        0.005277

I have read a paper that uses the CollectGcBiasMetrics on exome data and they magically produce a plot without bias. I say magically because they don't explain how they use CollectGcBiasMetrics.
Update I added the target coverage plot.

Update Am I missing something? Should I, for instance, use my target regions in instead of the reference genome when using CollectGcBiasMetrics? I will try it, but I am not sure this is the correct thing to do... It is hard to find any workflow for exome data with gatk and particularly CollectGcBiasMetrics, which specifies to use the reference genome.
From the comment @Maximilian Press

Do you have on-instrument diagnostics? Did you run FastQC on the raw reads? What do the other picard metrics look like? Did you manually inspect your alignments (strongly recommended)? Do you actually have coverage of the exome?

I will see if I can get some instrument diagnostics. I did look at the alignments, and as far as I can tell I don't see anything weird? Just that it is very deeply sequenced. I do see some intronic regions, with some reads, but only a few reads compared to the many on the targets, which I think is a consequence of the (too) deep sequencing. The fastqc of the raw reads didn't look troubling, and it did not give me any warnings. The picard metrics I thought did not look out of the ordinary. There is good target coverage. Yes, it is sequenced quite deep. Also, the lab tech did tell me they over-cycled because they wanted to be sure to have all the regions/targets.
Update So this bias might be due to too deep sequencing? And another important question, can I still reliably use this data? If so, I guess I can perform some corrections like suggested here
Any help/suggestions would be very welcome. Thank you in advance.

A bit of upvoting would be cool, its very good QC indeed and clearly presented. @supertech's answer makes sense to me. The GC% bias in living organisms could only represent bacteria IMO, I've seen this in soil bacteria so this has to be a lab artefact.

Thank you @M__. About the bacterial IMO (which I think means a contamination?). There is a separate run (another day a protocol slightly different) that shows the same bias. So either the lab tech is very sloppy, or the protocol was not correct (too many PCR cycles or something).

Thank you for incorporating the comments (and your answers to the comments) into your question. This is how comments should be used, and makes my job as a moderator much easier.

Am I missing something? Should I, for instance, use my target regions
instead of the reference genome?

Yes you should do this. Why do you think it's wrong doing that? It should be the first thing to check. And that GC bias might be effected by multiple factors. Most prominent one would be PCR over amplification. Most PCR polymerases give horrible GC bias when library is over cycled. There are many articles published on this.
I also suggest that you look at the number of reads in each bin. Peaks around 90% GC should not puzzle you so much. There is chance that you have only few reads in these bins out of millions of reads.

Just to formalise my comments.
@Supertech's is the right answer.
Background The one thing that is unusual is that the bias peaks at around 75%. The only time I've seen strong GC% bias is for particular groups of bacteria. 75% is about the maximum bias that can be achieved in a protein coding gene and these "extreme bacteria" hit that limit.
I better try and explain that a bit better. The primary selection pressure on a gene is the amino acid composition and the GC% is a second issue, i.e. selection pressure, for some organisms (e.g. certain bacteria, insects ..). Thus amino acids are not mutating, but the triplet codon degeneracy permits a bias in the GC% without alternating a given amino acids. In other words this GC% bias is occurring only for "neutral" mutations (clearly there must be a selection pressure, so they are not really neutral).
Summary and personal thoughts In summary the GC% limit is around 75% for a protein coding gene and beyond that it requires the amino acids to be changed. The peak in @Dandelion's data on that limit is weird because GC% artefacts from Taq errors are not subject to selection pressure: there should not be a peak.
"Contamination" I would avoid as an explanation because in my experience these groups of bacteria are unlikely to invade a lab environment. Sample mix up is a possible explanation. However, I would be cautious about assigning blame because in the chain of command within the wet-lab it could have occurred anywhere.
Official response To reiterate the formal and official answer would be @supertech's response and the 75% GC peak remains a mystery.

Please note The second plot has a distribution around GC content around 40% ... There appears a discrepancy between plots. Neither output appears compatible with the GC% of humans, but 40% GC with a reasonable proportion of that data on the 50% mark is better the 75% peak plot. Personally I would investigate the empirical basis for the differences.


text:
PLINK to VCF seems to to not convert correctly, what am I missing
I have Illumina Infinium Global SNP array data that I want to compare to same sample WES data in VCF format. To do that, I am trying to convert the plink files (.bed, .bim and .fam) into vcf.
I am using the following command:
plink2 --bfile gsa_prefix --fa {input.fa} --ref-from-fa --recode vcf id-paste=iid --out out_prefix

Somewhere I am missing something in the conversion or doing it wrong, it is not going properly. For some homozygous variants, the genotypes are switched, even though I think I tried to add the ref-from-fa to lock the reference to the REF alt, so I would expect this to be translated somehow in the conversion to vcf.
Here some examples (left before the || is the GSA, right after the || is WES):
sample 1:
chr3    10222796    rs3844283   C   G   .   .   ORIGINAL_CHROM=chr3;ORIGINAL_POS=10264480;ORIGINAL_STRAND=+;AC=2;AN=2   GT  1/1 || chr3 10222796    .   C   <*> 0   .   .   GT:GQ:MIN_DP:PL 0/0:50:17:0,51,509  106

sample2:
chr3    10222796    rs3844283   C   G   .   .   ORIGINAL_CHROM=chr3;ORIGINAL_POS=10264480;ORIGINAL_STRAND=+;AC=0;AN=2   GT  0/0 || chr3 10222796    .   C   G,<*>   70.9    PASS    .   GT:GQ:DP:AD:VAF:PL  1/1:67:104:0,104,0:1,0:70,70,0,990,990,990  99

So, this is not the case for all positions, some homozygous variants are the same genotype (i.e. 1/1 and 1/1). I know Illumina SNP array has this very elaborate way of annotating their "ALT"/"REF" alleles (or Allele A or B): https://www.illumina.com/documents/products/technotes/technote_topbot.pdf
But likely I am missing something because of this. I do not think this is technical sequencing wise (i.e. array just didn't pick up the SNP or the other way around). The pattern is too obvious over too many positions. And I just do not know now if a 0/0 means no variant/homozygeous ALT in terms of the reference genome or the other way around.
Note: It is a subsample of 22 samples from the plink files, it was subsampled by someone else.
Edit: could it be that my predecessor did not use --keep-allele-order when converting the array data to plink (with plink1.9)? The command that that I found in a random .txt file #ahlawdhelpme:
plink --file GSA_data --keep keep_list_22.txt --allow-no-sex --make-bed --out GSA_22_participants

I will look into this/try to redo this ^ and come with an update.
Edit2: the --keep-allele-order did not correct this.
Hope anyone has any insights. Thanks.


text:
true depth per gvcf location instead of min depth per gvcf block
I am using deepvariant to calls variants on my (WES) data. I am outputting a gvcf file.
I am wondering (I cannot find anything in the documentation) how it would be possible to get the true depth per site in the gvcf file and not only the min_dp per block. Currently I am converting the gvcf into vcf with bcftools for comparisons, but is there a way to calculate back what the true depth is for a site by adding some additional output to the gvcf file of some calculations?
Edit, the code I used:
/opt/deepvariant/bin/run_deepvariant --model_type={params.model} --ref={input.ref} --reads={input.bam} --output_vcf={output.vcf} --output_gvcf={output.gvcf} --num_shards {threads} --logging_dir={log}

bcftools convert --gvcf2vcf {input.gvcf} --fasta-ref {input.new_genome} --threads {threads} 2> {log.convert}

The file I get (in part):
chr1    69043   .   G   <*> 0   .   .   GT:GQ:MIN_DP:PL 0/0:12:4:0,12,119
chr1    69044   .   C   <*> 0   .   .   GT:GQ:MIN_DP:PL 0/0:12:4:0,12,119
chr1    69045   .   A   <*> 0   .   .   GT:GQ:MIN_DP:PL 0/0:12:4:0,12,119
chr1    69046   .   G   <*> 0   .   .   GT:GQ:MIN_DP:PL 0/0:12:4:0,12,119
chr1    69047   .   A   <*> 0   .   .   GT:GQ:MIN_DP:PL 0/0:12:4:0,12,119
chr1    21474900    .   A   G,<*>   64.3    PASS    .   GT:GQ:DP:AD:VAF:PL  0/1:63:91:49,41,0:0.450549,0:64,0,67,990,990,990

The non-variant sites have a MIN_DP that corresponds to the minimum depth in that block. So not per site. Compared to the variant that has a DP specific to that site.
Thanks in advance,

You would need to force an output for every site where you want depth information. if you just want depth, you can run samtools, but variant calling tools often take other things into account (base-quality, mapping, overlapping mates, etc) when calculating depth, so if you want that depth, with all the bells and whistles, you'll need to force an output on every base. I know this is possible with GATK, but I don't know about deep-variant.

I found that deepvariant does indeed not have this option. What I am looking for is apparently equivalent to the BP resolution option of GATK
-ERC BP_RESOLUTION.
The work around is to go into the bam files, and get the coverage per site. I am using bamtocov. But I might consider another one, since this one is extremely slow (even with multiple threads).
Edit: bamtocov took more than 6 hours and still was only at 1/4 through the bam file. I used samtools like Yossi Farjoun suggested.
samtools depth {input.bam} -b {input.targets} --threads {threads} --min-BQ {params.baseq} --min-MQ {params.mapq} -J -o {output}

It does have more options to control what is counted, like base quality, mapping quality, and if to include reads with deletions in depth computation.


text:
Variants from multiple tools: normalization before or after annotation with VEP?
We want to compare some variant calling tools and their calls on Whole Exome Sequencing data. We will have to normalize the variants that are called (.vcf format) before comparison using for instance vt normalize or hap.py. We are comparing with a gold standard too and use hap.py for this now to get some benchmark information.
Edit: with normalization I mean representing the variants in such a way, we are able to compare them: to ensure that the representation of a variant is “parsimonious” and “left-aligned.”
I was just wondering, and I cannot find any real guidelines (yet), when we annotate using VEP, do we have to normalize before or after annotation to keep the annotation correct as well. I know VEP has some normalization functionality too see here. So might we even not need an additional normalization tool before and just let VEP do it for us?
I am just wondering about the guidelines and what what others do.
Thanks in advance,

What kind of normalization are you going to do? VEP talks about normalizing variant represenations, variant names and left-aligning indels. However, all variant callers should be producing a vcf, so there is no need to normalize the names. Indels perhaps.  You might want to normalize other things like variant quality but we need to understand what you mean before we can help. Can you please [edit] your question and clarify?

If you're referring to left alignment + parsimonious representation, do it before any annotation. Most if not all annotation happens by comparing CHROM, POS and optionally REF and ALT and since vt normalize affects POS, you should do it before annotation.
I'm not sure how smart VEP is, but most annotators are simple lookup tools at least for third party databases, so they're not built to normalize while looking up.

Yes, that makes sense. And VEP has some form of normalization before annotation: here, but it seems to be more basic than a normalization tool like vt normalize.

Yes, vt normalize and bcftools norm both preserve a record of the pre-normalization entry so data loss is avoided. I'm not sure if VEP does that or even logs it.

Thinking a bit more about your comments ^, why would you not always normalize before annotation to be sure you would get the same results over different studies (hypothtically)? Would "the same variants" called with different callers, cause discrepencies in position and therefore annotation if you do not normalize them? Or am I missing something?

You're asking the right question. Without normalization, certain annotations could be missed especially if POS changes post normalization. Most callers should call a variant at the same location based on an alignment, but I'm not sure if that can be taken for granted. I always decompose and normalize all my VCFs pre annotation.

I guess for SNPs this might be (more) the case that they would have the same location between callers, INDELs are a different story likely

SNVs/MNVs that occur in a repeat region can confuse callers a little. One also needs to keep in mind that certain variants were named/annotated way back when and used to be called using the wrong coordinate for quite a while. I've seen this happen in old medical records, not sure how prevalent this is these days.


text:
Compare files with similar names using placeholders/wildcards in Snakemake
In Snakemake I want to use two input files for comparison. The one file has a sample name (e.g. 526846), with an experimental id (e.g. EF_S13' and 'MF_S5'). The file I want to compare it to has the same sample name only. I want to compare each sample with experiment ID (first files) to the file with corresponding sample name (second file). These files are in different directories. I thought I would be able to use two expand functions with two lists, as shown in the code below.
Note: these are example lists and the lists in reality is muchlarger.
Relevant code:

SAMPLE_SET = ['12345', '56789', '54321']
ID_SET = ['EF_S13', 'MF_S5', 'EF_S11', 'MF_S3', 'EF_S10', 'MF_S2']

rule all:
    input:
        array = expand("results/{sample}_{samp_id}_array/{sample}_{samp_id}.summary.csv", zip, sample = SAMPLE_SET, samp_id = ID_SET)

rule happy_rtg_array:
    """ Normalizes and compares vcf files with array samples"""
    input:
        query_vcf = expand("WESP_run/results/variant_calling/{sample}_{samp_id}.vcf", zip, sample = SAMPLE_SET, samp_id = ID_SET),
        truth_vcf = expand(os.path.join(config["truth_dir_vcf_gsa"],"{sample}.vcf.gz"), sample = SAMPLE_SET),
    output:
        "WESPValidation/results/{sample}_{samp_id}_array/{sample}_{samp_id}.summary.csv"
[more code]

However, when I look at the output log file with --dry-run, I get as input for the job a huge list with all the inputs in one go, rather then per file/comparison, and one output file (see relevant output), and that times the amount of the (right) combinations of the file name and id (counted from the output).
Also the wildcards are not what I thought they should be. I might be confused with the expand function, as I thought how I used them now should "lock in" the wildcards and become placeholders filled in by the lists.
Relevant output:
[Tue Oct  3 14:31:20 2023]
rule happy_rtg_array:
    input: 12345_EF_S13.vcf, 12345_MF_S5.vcf, 56789_EF_S11.vcf, 56789_MF_S3.vcf, 54321_EF_S10.vcf, 54321_MF_S2.vcf, 12345.vcf, 56789.vcf, 54321.vcf
    output: results/12345_EF_S13_array/12345_EF_S13.summary.csv
    log: happy_rtg_array/12345_EF_S13.log
    jobid: 6
    reason: Missing output files: WESPValidation/results/12345_EF_S13_array/12345_EF_S13.summary.csv
    wildcards: sample=12345_EF, samp_id=S13
    threads: 10
    resources: mem_mb=10000, disk_mb=7855, tmpdir=/tmp, threads=1, runtime_min=120

I doubt this is how it should be, and wonder what I am doing wrong. Is there a better/correct way to do this.
Thank you in advance

After some trial and error, I fixed it with the expand function only in the rule all and removing the expand from the input in the rule happy_rtg_array (thanks to @Pallie), and by changing the output file paths:
SAMPLE_SET = ['12345', '56789', '54321']
ID_SET = ['EF_S13', 'MF_S5', 'EF_S11', 'MF_S3', 'EF_S10', 'MF_S2']

rule all:
    input:
        array = expand("results/{sample}/{samp_id}/summary.csv", zip, sample = SAMPLE_SET, samp_id = ID_SET)

rule happy_rtg_array:
    """ Normalizes and compares vcf files with array samples"""
    input:
        query_vcf = "_".join(["WESP_run/results/variant_calling/{sample}","{samp_id}.vcf"]),
        truth_vcf = os.path.join(config["truth_dir_vcf_gsa"],"{sample}.vcf.gz"),
    output:
        "WESPValidation/results/{sample}/{samp_id}/summary.csv"
[more code]

Basically, I replaced all the underscores with slashes. For some reason these underscores are messing up the wildcards when trying to connect them with it. Maybe someone more knowledgeable knows why this is.
[Tue Oct  5 12:24:10 2023]
rule happy_rtg_array:
    input: 12345_EF_S13.vcf
    output: results/12345/EF_S13/summary.csv
    log: happy_rtg_array/12345_EF_S13.log
    jobid: 6
    reason: Missing output files: WESPValidation/results/12345/EF_S13/summary.csv
    wildcards: sample=12345, samp_id=EF_S13
    threads: 10
    resources: mem_mb=10000, disk_mb=7855, tmpdir=/tmp, threads=1, runtime_min=120

For the reason why, check the part about ambiguous wildcards: https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#wildcards. I missed the amibguity of your wildcards when I wrote my answer :)

Dank je @Pallie :-)

Using expand in your rule's input:  causes snakemake to aggregate all vcfs into one list and leads to your undesired behaviour. By simply removing the expand calls from the input: you should get your desired result.
https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#the-expand-function

This is only partially the answer, I tried this before indeed, but still get the following wildcards "wildcards: sample=12345_EF, samp_id=S13" and they should be: "wildcards: sample=12345, samp_id=EF_S13", to get the right file name for the input "truth_vcf"


text:
VEP annotation INFO field Ensembl IDs and locations
I have a vcf file that I annoteted with VEP, for human data. I have run VEP to annotate my files with some additional parameters (as shown below in the ##VEP-command-line).
However, my output is rather strange (mainly the INFO column).
##VEP="v108" time="2023-04-27 15:13:08" cache="workflow/resources/variants/cache_vep/homo_sapiens/108_GRCh38" ensembl-funcgen=108.56bb136 ensembl-variation=108.a885ada ensembl-io=108.58d13c1 ensembl=108.d8a9c80 1000genomes="phase3" COSMIC="96" ClinVar="202205" HGMD-PUBLIC="20204" assembly="GRCh38.p13" dbSNP="154" gencode="GENCODE 42" genebuild="2014-07" gnomADe="r2.1.1" gnomADg="v3.1.2" polyphen="2.2.2" regbuild="1.0" sift="sift5.2.2"
##INFO=<ID=CSQ,Number=.,Type=String,Description="Consequence annotations from Ensembl VEP. Format: Allele|Consequence|IMPACT|SYMBOL|Gene|Feature_type|Feature|BIOTYPE|EXON|INTRON|HGVSc|HGVSp|cDNA_position|CDS_position|Protein_position|Amino_acids|Codons|Existing_variation|DISTANCE|STRAND|FLAGS|SYMBOL_SOURCE|HGNC_ID|GENE_PHENO|SIFT|PolyPhen|AF|AFR_AF|AMR_AF|EAS_AF|EUR_AF|SAS_AF|gnomADe_AF|gnomADe_AFR_AF|gnomADe_AMR_AF|gnomADe_ASJ_AF|gnomADe_EAS_AF|gnomADe_FIN_AF|gnomADe_NFE_AF|gnomADe_OTH_AF|gnomADe_SAS_AF|CLIN_SIG|SOMATIC|PHENO">
##VEP-command-line='vep --dir_cache workflow/resources/variants/cache_vep --vcf -i WESP_run/results/variant_calling/VP_DBS_25_7_S7.vcf -o WESP_run/results/variant_annotation/VP_DBS_25_7_S7.annotated.vcf.gz --af --af_1kg --af_gnomade --sift b --polyphen b --gene_phenotype --fork 24 --offline --compress_output gzip'
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  VP_DBS_25_7_S7
chr1    13273   .   G   C   25.3    PASS    CSQ=C|non_coding_transcript_exon_variant|MODIFIER|DDX11L1|ENSG00000223972|Transcript|ENST00000450305|transcribed_unprocessed_pseudogene|5/6||||313|||||rs531730856||1||HGNC|HGNC:37102||||0.0950|0.0204|0.1455|0.0625|0.1471|0.1401||||||||||||,C|non_coding_transcript_exon_variant|MODIFIER|DDX11L2|ENSG00000290825|Transcript|ENST00000456328|lncRNA|3/3||||521|||||rs531730856||1||EntrezGene|||||0.0950|0.0204|0.1455|0.0625|0.1471|0.1401||||||||||||,C|downstream_gene_variant|MODIFIER|WASH7P|ENSG00000227232|Transcript|ENST00000488147|unprocessed_pseudogene||||||||||rs531730856|1131|-1||HGNC|HGNC:38034||||0.0950|0.0204|0.1455|0.0625|0.1471|0.1401||||||||||||,C|downstream_gene_variant|MODIFIER|MIR6859-1|ENSG00000278267|Transcript|ENST00000619216|miRNA||||||||||rs531730856|4096|-1||HGNC|HGNC:50039||||0.0950|0.0204|0.1455|0.0625|0.1471|0.1401|||||||||||| GT:GQ:DP:AD:VAF:PL  0/1:23:55:28,27:0.490909:25,0,27

For instance, at the first location/variant I look at the Ensembl gene id and see multiple genes on the same locations and strand. This is very strange to me, as I would understand different transcripts, but different genes on the same strand? Also, apart from that, the corresponding locations/regions of genes on the Ensembl website do not agree with the location in the vcf file.
For instance, the location for this SNP is 13273, but the annotations are (corresponding ensembl archive version 108):

Agrees: http://oct2022.archive.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=ENSG00000223972;r=1:12010-13670;t=ENST00000450305
Agrees: http://oct2022.archive.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=ENSG00000290825;r=1:11869-14409;t=ENST00000456328
Does not agree: http://oct2022.archive.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=ENSG00000227232;r=1:14404-29570;t=ENST00000488147
Does not agree: http://oct2022.archive.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=ENSG00000278267;r=1:17369-17436;t=ENST00000619216

What am I missing here?
Kind regards,
Eva

You need to set your expectations clearly instead of just saying "Does not agree". I'm taking a guess here - your statement is that the 2/4 of the gene co-ordinates do not overlap the variant location.
This is because you're viewing just part of the picture - the co-ordinates - and not seeing another important part, which is the kind of effect. For the latter two genes, VEP clearly says the effect is downstream_gene_variant, which Sequence Ontology defines as "A sequence variant located 3' of a gene."
Putting those two pieces together, the entire picture makes sense. Variants affect different transcripts (of possibly different genes) in different manners. Looking at just 1 of those 2 components will of course make things seem odd.

You are right, it is just overlapping genes. And VEP annotates these, while other variant annotation tools like ANNOVAR do not (just outputs the most deleterious) ref.

That’s not true. Both tools can annotate all or most deleterious effects depending on the parameters you use.

Just saw the reference post you’d added - that post is wrong. See ANNOVAR documentation- it says “ If the users want to have all functional consequences printed out (rather than just the most important one defined by the precedence above), the --separate argument should be used. In this case, several output lines may be present for each variant, representing several possible functional consequences.”

I guess the difference would be the default behavior maybe, but you are right. Sorry for the confusion and thank you for your patience.


text:
Stringtie/DEprep.py gene/transcript IDs are wrongly formatted
Hi my RNAseq workflow is ending up with wrongly formatted gene IDs (and separately transcript _IDs) after a Hisat2 ->samtools sort -> stringtie -> DEprep.py workflow.
DEprep.py outputs a gene_counts.csv which includes this...
gene:ENSRNOG00000000017         47
ENSRNOG00000000021              220

Spelling this out "gene:" should not be in the gene_ID.  First thought of course was to just delete the "gene:" string with search and replace.
However, sometimes, the same gene appears as both indicating the error could be really damaging:
Gene_ID                   Counts
ENSRNOG00000004937        5234
gene:ENSRNOG00000004937   0

Digging around I have found that the issue is present in the initial stringtie /ballgown/ GTF output prior to running DEprep.py
chr1    ensembl exon    la la   gene_id
 "gene:ENSRNOG00000046319" la la la

Curiously my workflow is copied from my own local Galaxy installation workflow that uses the same steps together with the SAME fastq's the same, Hisat2 index, the same Guide.GFF, but there is no issue there.
I have spent a day googling this and no success yet, lots of peep have missing ids. One other person had this symptom, but sadly no resolution.  Any suggestions gratefully received!!
Code details below and an attached folder in the comments with dummy data and r script that reproduces this error.
HiSat2:

hisat2 --dta --add-chrname --time --summary-file YH2minihisat.txt
--rna-strandness FR  --pen-cansplice 0 --pen-noncansplice 12 --pen-canintronlen G,-8.0,1.0 --pen-noncanintronlen G,-8.0,1.0 --min-intronlen 20 --max-intronlen 500000 --threads 12 -x rn6/genome -q -1 YH2mini_R1.fastq.gz -2 YH2mini_R2.fastq.gz |  samtools sort -m20G -l 0 -T /Users/me/Documents/TEMP/b -O bam | samtools view -O bam -@ 12 -o YH2mini.bam

Stringtie:

stringtie -fr -p 12 -G rat6.gtf -e -b ballgown/YH2mini -f 0.15 -m 200
-a 10 -j 1 -c 2 -g 50 -M 0.95 -o ballgown/YH2mini/YH2mini.gtf YH2mini.bam

PrepDE:

prepDE.py -i ballgown -g YH2mini_gene_counts.csv -t
YH2mini_transcript_counts.csv -l 100

So peep can see I'm not lazy, I have a work around where I can run this, extra prior to prepDE.py, but seems very wrong!!
if(CORRECTGTF){

        command <- paste0("sed 's/gene\\://g' ballgown/",PREFIX,"/",PREFIX,".gtf > temp.txt")
        system(command)
        command <- paste0("sed 's/transcript\\://g' temp.txt > ballgown/",PREFIX,"/",PREFIX,".gtf")
        system(command)
        command <- paste0("sed 's/gene\\://g' ballgown/",PREFIX,"/","t_data.ctab > t_datafix.ctab.tmp")
        system(command)
        command <- paste0("sed 's/transcript\\://g ' t_datafix.ctab.tmp > ballgown/",PREFIX,"/","t_data.ctab")
        system(command)
    }

The gene:ENSRNOG00000046319 in the GTF file is pretty standard. I don't think we can help without access to the actual data. Can you give us some examples? for instance, you seem to be searching through something for the gene names. What are you searching through? Presumably the output of your pipeline, but what kind of output is that? Then, you say that you sometimes see both but you show a count of 0 for gene:ENSRNOG00000004937. My guess is that you are getting the gene names from different fields, but I would need to understand what you're doing in more detail to be sure.

@ Terdon Thanks, I will just need to investigate now how to share data and upload anything needed (after a meeting now!).

We don't need real data, per se, just some examples to understand what you mean.

...OK, sorry... don't really know how to give examples beyond the one above and so the actual code and data are here... https://www.dropbox.com/sh/lqk0usgt8n2suvx/AADqZDNpKKPKFGWZQ9sn9ix6a?dl=0


text:
How to adjust #columns shown when PyRange is printed?
In Pandas, we have pandas.options.display to adjust how a DataFrame should be displayed when printed.
I wonder if PyRanges has a similar set of options?
For example, if I have a PyRange object like this:
import pyranges as pr
data = pr.data.ucsc_bed()

I can do print(data) and get:
+--------------+-----------+-----------+------------+------------+-------+
| Chromosome   | Start     | End       | Feature    | gene_id    | +4    |
| (category)   | (int32)   | (int32)   | (object)   | (object)   | ...   |
|--------------+-----------+-----------+------------+------------+-------|
| chr1         | 12776117  | 12788726  | gene       | AADACL3    | ...   |
| chr1         | 169075927 | 169101957 | gene       | ATP1B1     | ...   |
| chr1         | 6845383   | 7829766   | gene       | CAMTA1     | ...   |
| chr1         | 20915589  | 20945396  | gene       | CDA        | ...   |
| ...          | ...       | ...       | ...        | ...        | ...   |
| chrX         | 152661096 | 152663330 | exon       | PNMA6E     | ...   |
| chrX         | 152661096 | 152666808 | transcript | PNMA6E     | ...   |
| chrX         | 152664164 | 152664378 | exon       | PNMA6E     | ...   |
| chrX         | 152666701 | 152666808 | exon       | PNMA6E     | ...   |
+--------------+-----------+-----------+------------+------------+-------+
Stranded PyRanges object has 5,519 rows and 9 columns from 30 chromosomes.
For printing, the PyRanges was sorted on Chromosome and Strand.
4 hidden columns: transcript_id, Strand, exon_number, transcript_name

How can I increase the number of columns shown?
Also appreciated if you refer me to a documentation that explains other related display options.
Update: For future reference, @terdon remark encouraged me to check if such an adjustment is even considered in the PyRanges. I found the line that is responsible for this representation, and indeed the COLUMN environment variable is the only defining parameter. In conclusion, the adjustment that I am looking for is simply not implemented in PyRanges.

Are you running in an interactive python shell? On what operating system? What happens if you make the terminal window larger? Or if you print to a file? Can you provide a minimal code example that we can run which reproduces the problem?

Thanks @terdon for your reply. I already provided the minimal code example. Just run the code I provided and you get this picture.

The only code you have provided is data = pr.data.ucsc_bed(). We would need an example file, and the commands to read it into the script etc. Please see https://stackoverflow.com/help/minimal-reproducible-example. Also, please answer my other questions too.

@terdon, please note that ucsc_bed() produces a dummy data. I now added import pyranges as pr. Just copy-paste this code and you get the output. And regarding your questions: I am answering :). Sorry, I pressed Enter and the comment got posted. I tested this in the terminal and also jupyter notebook, but the output is independent of the environment. It also stays the same irrespective of the size of window.

Ah! Thank you, I hadn't realized that. And thanks for the answers! Does it print the whole thing if printed to a file, at least?

Of course. You can always get the Pandas.DataFrame using data.df, and then display it as you wish. Alternatively, you can export the bed file using data.to_bed('./test.bed'). In all cases, all data is exported. Note that my question is related to display of the data in the terminal. Not exporting issues. Please consider taking a look at pandas.options.display to understand what I mean by difference between exporting the PyRanges objects, and displaying PyRanges objects.

Let us continue this discussion in chat.

When trying this on my machine, if I resize the terminal window before running print(data), I get the full output:

You can also force this, even if your terminal is actually smaller, by setting the COLUMNS variable to something larger than the actual size when layching python or ipython. This is a special bash variable that holds the current size of the terminal. From man bash:
   COLUMNS
          Used  by  the select compound command to determine the terminal
          width when printing selection lists.  Automatically set if  the
          checkwinsize  option is enabled or in an interactive shell upon
          receipt of a SIGWINCH.

So if you set it manually, you can get the entire data frame printed, BUT you don't really want to do that since it won't scroll, it just wraps the text. For example, try:
COLUMNS=300 python

That will set the size to 300, which results in this:

Very good @terdon. Thanks! I tested this with Jupyter notebook and even when adjusting the window size following this the additional columns are still kept hidden. Most likely, I should change another parameter in Jupyter. But indeed, I could confirm that in Terminal if I adjust the window size (even inside a screen session), the columns do show as long as the window size is large enough. Your approach using COLUMN variable is very nice!

You're welcome. If Jupyter gives you a shell, you can try setting the variable in the shell. Alternatively, you can actually do it in your python code with os.environ['COLUMNS']="300" and then run print(data) and it might work as expected (confirmed in an ipython session, but I don't have a jupyter notebook to test on).

I confirm that setting the environment variable using os.environ['COLUMNS']="300" works in a Jupyter notebook too! Thank you very much @terdon.


text:
extracting raw gene counts based on the cluster generated by seurat
I'm new to single cell sequencing analysis. I have made different clusters using seurat. Now I would like to extract the raw gene counts based on the cluster generated by seurat.
i have provided my R code used in the analysis.
library(dplyr)
library(Seurat)
library(patchwork)
ftp="https://ftp.ncbi.nlm.nih.gov/geo/series/GSE201nnn/GSE201254/suppl/GSE201254_count_mat_Control-NoCrush.csv.gz"
counts=data.table::fread(ftp)
C_nocrush=counts %>%
select(contains(c("V1","CtCReRGCS")))
rownames(C_nocrush)=C_nocrush$V1
C_nocrush$V1=NULL
CNC=CreateSeuratObject(counts = C_nocrush,project = "He RGC Data", min.cells = 3, min.features = 200)
CNC[["percent.mt"]] = PercentageFeatureSet(CNC, pattern = "^MT-")
VlnPlot(CNC, features = c("nFeature_RNA", "nCount_RNA","percent.mt"), ncol = 3)
plot1 <- FeatureScatter(CNC, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(CNC, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1+plot2
CNC = subset(CNC, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
CNC = NormalizeData(CNC, normalization.method = "LogNormalize", scale.factor = 10000)
CNC = FindVariableFeatures(CNC, selection.method = "vst", nfeatures = 2000)
top10 <- head(VariableFeatures(CNC), 10)
plot1 <- VariableFeaturePlot(CNC)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1+plot2
all.genes <- rownames(CNC)
CNC <- ScaleData(CNC, features = all.genes)
CNC <- RunPCA(CNC, features = VariableFeatures(object = CNC))
VizDimLoadings(CNC, dims = 1:2, reduction = "pca")
DimPlot(CNC, reduction = "pca")
DimHeatmap(CNC, dims = 1:15, cells = 500, balanced = TRUE)
CNC <- JackStraw(CNC, num.replicate = 100)
CNC <- ScoreJackStraw(CNC, dims = 1:20)
JackStrawPlot(CNC, dims = 1:15)
ElbowPlot(CNC)
CNC <- FindNeighbors(CNC, dims = 1:10)
CNC <- FindClusters(CNC, resolution = 0.5)
head(Idents(CNC), 5)
CNC <- RunUMAP(CNC, dims = 1:10)
DimPlot(CNC, reduction = "umap")
saveRDS(CNC, file = "Desktop/Single_cell_RNA/CNC.rds")
cluster2.markers <- FindMarkers(CNC, ident.1 = 2, min.pct = 0.25)
head(cluster2.markers, n = 5)
new.cluster.ids=c("T5/F RGCs", "T5/F RGCs","T-RGCs","Other RGcs and cartpt RGCs","T5/F RGCs","IP-RGCs","Other RGcs and cartpt RGCs","N-RGCs","Alpha-RGCs" )
names(new.cluster.ids)=levels(CNC)
CNC=RenameIdents(CNC,new.cluster.ids)
 DimPlot(CNC,reduction = "umap", label = T, pt.size = 0.5)+NoLegend()

Tried using this command to extract the raw counts based on the cluster but is not working for me.
meta.data.cluster=unique(x=CNC@meta.data$active.ident)
meta.data.cluster=unique(Idents(CNC))
for(group in meta.data.cluster) {
    group.cells <- WhichCells(object = CNC, idents = group)
    data_to_write_out <- as.data.frame(GetAssayData(CNC,slot = 'counts')[,group.cells])
    write.csv(data_to_write_out, row.names = TRUE, file = "1_cluster_outfile.csv")
}

Hi, thank you for supplying detailed code. Could you kindly explain the precise question? I assume the code is throwing an error or not working how it was anticipated. Describing this would be great.

I have updated my script. I want to extract the raw counts of the genes using the selected cluster.

How do you know that the count extraction code you wrote is not working?


text:
Is it common to get different number of SNVs+Indels across samples from vcf files generated using GATK and DRAGEN (counts are higher for GATK)?
We have 2 vcf (Whole Exome Sequencing (WES) data; germline samples) files (e.g., vcf_1 and vcf_2). vcf_1 was generated (Ref. genome: hg38) using the GATK pipeline for 250 children and their parents and vcf_2 was generated (Ref. genome: hg38) using the DRAGEN pipeline for another 50 children and their parents (Samples are different). These are raw vcf files and I have not applied any quality-related false positive filters to them.
Now, when we compare the sample-wise SNV counts from vcf_1 and vcf_2, there is a huge difference. vcf_2 contains on average 4-6 times lesser number of SNVs per sample than vcf_1.
N.B:
I found this article (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4418901/), where it has been reported that the mean numbers of single-nucleotide variants (SNVs) and small insertions/deletions (indels) detected per sample was around 98k, respectively, for WES.  For our case, in the case of GATK-generated vcf samples, this varies from 130k to 188k, and for DRAGEN-generated vcf samples, this varies from 30k to 40k.
I understand both pipelines (GATK and DRAGEN) implement different methods to generate vcf files but is it common/normal to examine that much difference in SNV counts? I would appreciate it if someone can give some explanation (or provide some materials) about it.

We need a lot more detail to be able to help. What algorithm was used? What version of the various tools? Are these germline samples or somatic? What species? What kind of sample (WGS, WES, panels)? What filtering options were applied? Please [edit] your question and add more detail.

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

I have updated the question with a few details (WES data, human sample, Ref. genome used: hg38, etc.) as you have mentioned. As soon as I will get some more info about the commands used to generate the vcf files or the version of GATK and DRAGEN, I will update the question further. Many thanks in advance.

I suspect what's happening here is an unfortunate side-effect of making variant-finding algorithms memory-efficient, so that they can work on extremely large populations. By doing the processor-intensive variant calling separately on each individual, calculations can be easily parallelised to run quickly.
For the numbers you have (i.e. 300 children and their parents), there's an alternative approach for generating variant calling that uses all mapped files as input, rather than only one at a time, i.e. instead of running:
bcftools mpileup [options] in1.bam
bcftools mpileup [options] in2.bam
bcftools mpileup [options] ...

The all-at-once approach has a single mpileup call:
bcftools mpileup [options] in1.bam in2.bam ...

This will call multiple samples at the same time, so that the output will have identical variant lines for all individuals. It also has the added benefit of getting more information on variation, allowing more rare variants to be identified.
I'm not sure which approach is being used for DRAGEN vs GATK, but I do know that at some stage in the past the GATK standard practice was calling one sample at a time. Curiously, if this were still the case, I'd expect the results to be the other way round, as the one-at-a-time approach tends to be less sensitive.
I guess it's also possible that the pipelines have different sensitivity / specificity tweaks, that could also explain the differences in variant counts.


text:
Get genome coordinates from UniProt entry?
I'm trying to get genome coordinate information programmatically for multiple genes.
As an example, I'd like to get the annotation (list of exons and their chromosomes, starts, and ends) for this UniProt entry: https://www.uniprot.org/uniprotkb/B6TU39/entry.
There is a corresponding GenBank entry, which gives the CDS (https://www.ncbi.nlm.nih.gov/nuccore/EU968504) but I can't see anywhere how to extract the genome coordinates. How can I do this?
Thanks!

You have asked for a 'programmatic' solution; what programming language?

I'll be using Python, but anything using the UniProt/GenBank/etc. REST APIs would work.

In the first instance it'd be useful just to know how to do it using the website GUIs, as that's likely to give enough clues

Talking of web GUIs, you can use Ensembl plants. Put the uniprot accession in the search box and you will find the gene page.

You don't want the EU968504, that's the mRNA not the genome sequence. You want https://www.ncbi.nlm.nih.gov/nuccore/CM007647, also linked to from the UniProt page, just under the link to EU968504:

You also need the gene name, ZEAMMB73_Zm00001d029604 in this case, available further up on the page in the "Names & Taxonomy" section:

Next, we can download the GenBank file for the relevant chromosomal contig  using NCBI's entrez API:
curl "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=CM007647&rettype=gbwithparts&retmode=text" > CM007647

This will take a few minutes. Once that is done, you can parse the file and extract the start and end coordinates using something like this:
#!/usr/bin/env python3
from Bio import SeqIO
import sys

gb_record=SeqIO.read(open(sys.argv[1],"r"), "genbank")

for feature in gb_record.features:
    if feature.type == "gene" and sys.argv[2] in feature.qualifiers['locus_tag']:
        if feature.location.strand == 1:
            strand = "+"
        elif feature.location.strand == -1:
            strand = "-"
        else:
            print("Unknown strand: ", feature.location.strand)
            sys.exit(1)

        print("%s\t%s\t%s\t%s\t%s\t%s\t%s\t.\t." % (gb_record.id, "GenBank",
                                                    feature.type, feature.location.start,
                                                    feature.location.end, ".", strand))

You pass the GenBank file and gene name as arguments, and this is the result:
$ foo.py CM007647  ZEAMMB73_Zm00001d029604
CM007647.1  GenBank gene    78770080    78773320    .   +   .   .

Thanks so much for the detailed answer!

@MichaelDunne CM007647 comes from GRAMENE-4.0, not the latest NAM-5.0 maize genome. On NAM-5.0, this gene is located at 1:77,982,020-77,986,746.

You could align back the peptide sequence to the genome, respecting splicing, for example using miniprot from Heng Li.

As it’s currently written, your answer is unclear. Please [edit] to add additional details that will help others understand how this addresses the question asked. You can find more information on how to write good answers in the help center.

Hi and welcome to the site! We like answers to be a bit more comprehensive here. Could you expand your answer to explain how to use miniprot to do what the OP asks for? Also, you would need to explain how to go from a protein accession to its fasta sequence to then use that to align to the genome, as well as how to figure out which genome to align to and how to differentiate between homologs.


text:
How to handle output from a python script that gives different outputs depending upon the condition - Nextflow
I have created three Nextflow processes: One checks the strandedness of the fastq file and outputs it to a .txt file. Second one executes a python script with different commands depending upon the condition. And the third one aligns them.
/*
 * the `CHECK_STRANDEDNESS` process takes SAM file
 * and a BED file
 * and outputs a a txt file containing information about strandedness.
 */

process CHECK_STRANDEDNESS {

    input:
    path "${pair_id}.sam"
    path "${gtf.baseName}.bed"

    output:
    path "${pair_id}_strand_check.txt"

    script:
    """
    infer_experiment.py -i "${pair_id}.sam" -r "${gtf.baseName}.bed" > "${pair_id}_strand_check.txt"
    """
}

/*
 * the `FINDTAIL` process takes "${pair_id}_strand_check.txt" file
 * and executes commands according to the reads strandedness.
 * The output generated is FASTA files.
 */

process FINDTAIL {

    input:
    path "${pair_id}_strand_check.txt"
    path tuple val(pair_id), path(reads)

    output:

    script:
    """
    #!/usr/bin/env python3

    import subprocess
    import pandas as pd
    import os

    result = pd.read_csv("${pair_id}_strand_check.txt", sep="\r\n", header=None, engine='python')

    failed = float(result.iloc[1,0].replace('Fraction of reads failed to determine: ', ''))
    fwd = float(result.iloc[2,0].replace('Fraction of reads explained by "1++,1--,2+-,2-+": ', ''))
    rev = float(result.iloc[3,0].replace('Fraction of reads explained by "1+-,1-+,2++,2--": ', ''))
    fwd_percent = fwd/(fwd+rev)
    rev_percent = rev/(fwd+rev)

    if float(result.iloc[1,0].replace('Fraction of reads failed to determine: ', '')) > 0.50:
        cmd1 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_1.fasta
        cmd2 = "findtail_v1.01 --input_file ${reads[1]}  --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_2.fasta
        print(cmd1)
        subprocess.call(cmd1, shell=True)
        print(cmd2)
        subprocess.call(cmd2, shell=True)

        if fwd_percent > 0.9:
            #Over 90% of reads explained by "1++,1--,2+-,2-+
            #Data is likely FR/fr-secondstrand
            cmd1 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > tail_${pair_id}_1.fasta
            cmd2 = "findtail_v1.01 --input_file ${reads[1]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > tail_${pair_id}_2.fasta
            print(cmd1)
            subprocess.call(cmd1, shell=True)
            print(cmd2)
            subprocess.call(cmd2, shell=True)

        elif rev_percent > 0.9:
            # Over 90% of reads explained by "1+-,1-+,2++,2--
            # Data is likely RF/fr-firststrand
            cmd1 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_1.fasta
            cmd2 = "findtail_v1.01 --input_file ${reads[1]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_2.fasta
            print(cmd1)
            subprocess.call(cmd1, shell=True)
            print(cmd2)
            subprocess.call(cmd2, shell=True)

        elif max(fwd_percent, rev_percent) < 0.6:
            #Under 60% of reads explained by one direction
            #Data is likely unstranded
            cmd1 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_1.sl.fasta
            cmd2 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > tail_${pair_id}_1.pr.fasta
            cmd3 = "findtail_v1.01 --input_file ${reads[1]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_2.sl.fasta
            cmd4 = "findtail_v1.01 --input_file ${reads[1]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > tail_${pair_id}_2.pr.fasta
            print(cmd1)
            subprocess.call(cmd1, shell=True)
            print(cmd2)
            subprocess.call(cmd2, shell=True)
            print(cmd3)
            subprocess.call(cmd3, shell=True)
            print(cmd4)
            subprocess.call(cmd4, shell=True)

        else:
            #if strand couldnt be detected
            # we assume it is first-stranded
            # since most illumina reads are
            cmd1 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_1.fasta
            cmd2 = "findtail_v1.01 --input_file ${reads[1]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_2.fasta
            print(cmd1)
            subprocess.call(cmd1, shell=True)
            print(cmd2)
            subprocess.call(cmd2, shell=True)
        """
}

/*
 * define the `ALIGN` process that outputs unmapped files in FASTA format
 * It will take the FASTA files creates by the `FINDTAIL` process
 * and try to align them INDIVIDUALLY.
 */

process ALIGN {

    input:
    path index
    path "tail_${pair_id}*.fasta"

    output:
    path "unmapped_tail_${pair_id}*.fasta"

    script:
    """
    bowtie2 -x $index -f "tail_${pair_id}*.fasta" --un "unmapped_tail_${pair_id}*.fasta"
    """
}

After creating a "${pair_id}_strand_check.txt" file, I input it to the FINDTAIL process that finds reads with polyadenylated ends. It takes fastq files individually, so I have to run it individually for each file.
If fastq file pair is first stranded or second stranded, I executed findtail twice,once for R1 and once for R2, which creates two fasta files. But if it is unstranded, I have to run findtail four times, twice for each file. Thus, depending on the paired-end reads' strandedness, I have different number of inputs. How do I define the output of this process?
In the next process I align these files individually to get unmapped reads.
This issue has been eating my brain away for 5 days. So please let me know if something is unclear.
Best wishes!

From the (now removed) edits, the updated path glob pattern could just be: path("tail_${pair_id}_{1,2}.{pr,sl}.fasta"). If it's still not working, check to see what files are actually being created. They might be being created in another directory or not at all.

I ran it the same thing again and it worked somehow. So  gathered it was fault on my end and remove the edit.

You should be able to just use a glob pattern that accounts for all possible output files, regardless of the specified conditions. As long as a single file matches the glob pattern when the process completes (successfully), Nextflow won't complain. I think this should suffice, but please let me know in the comments if it doesn't:
output:
tuple val(pair_id), path("tail_${pair_id}_{1,2}{,.pr,.sl}.fasta")

Thank you so much! Merry Christmas and Happy Holidays!


text:
How to manually trim FASTA file sequences with the information provided in the header and store it into a new FASTA file? - Python
I have a FASTA file that has reads like this:

>SRR5655563.745 745 length=126 (type=T,start=2,end=12,length=11,identity=81.8182%) CttttatttgttGTGTGTGAGGTTTGATTTGATGGGAAAATATCTTGAATCTGCGGCGAGGTTGGAAGAGCTATCGCGGATTGTGTCATCTGCTGCGAAGCCCAATAGGTCAAAGGGAACGCTACC

As you can see, there is information about the start "start=2" and end "end=12" of a sequence within the header. I would like to slice the sequence like [start:end] and keep the rest of it. e.g
The part I would like to trim:

Ctttggtttcctttt

And after trimming, I would like to keep the rest of the read:

GTGTGTGAGGTTTGATTTGATGGGAAAATATCTTGAATCTGCGGCGAGGTTGGAAGAGCTATCGCGGATTGTGTCATCTGCTGCGAAGCCCAATAGGTCAAAGGGAACGCTACC

So the new FASTA file will have reads like this:

>SRR5655563.745 745 length=126 (type=T,start=2,end=12,length=11,identity=81.8182%)
GTGTGTGAGGTTTGATTTGATGGGAAAATATCTTGAATCTGCGGCGAGGTTGGAAGAGCTATCGCGGATTGTGTCATCTGCTGCGAAGCCCAATAGGTCAAAGGGAACGCTACC

I have written a python script that takes a FASTA file as an input and outputs a trimmed FASTA file, like this:
tailtrim.py --fasta input.fasta --output output.tail.trimmed.fasta

The code does trim the reads where I want it to.
import numpy as np
import argparse
import re

def main():
    parser = argparse.ArgumentParser()

    parser.add_argument("--fasta", type=str, required=True, help="input fasta file")
    parser.add_argument("--output", type=str, required=True, help="output tail.trimmed.fasta file")

    args = parser.parse_args()

    out_fasta = open(args.tail.trimmed.fasta, 'w')

    for line in open(args.fasta):
        header = line.rstrip()
        read = next(args.fasta).rstrip()

    #Extract header into a list
        header_lst = []
        header_lst= header.split(" ")

        tmp_start_lst = []
        tmp_end_lst = []
        tmp_lst = []

        for i in header_lst:
            if i.startswith("(") and i.endswith(")"):
                tmp_lst.append(i.split(","))
        #print(tmp_lst)

    # Extract the start and end numbers in a list
        for i in tmp_lst:
            for str in i:
                if str.startswith("start="):
                    tmp_start_lst.append(str)
                if str.startswith("end="):
                    tmp_end_lst.append(str)

        pattern = re.compile(r'\d+$')

        start_lst = []
        for str in tmp_start_lst:
            match = re.search(pattern, str)
            start_lst.append(int(match.group()))

        end_lst = []
        for str in tmp_end_lst:
            match = re.search(pattern, str)
            end_lst.append(int(match.group()))

        #Extract reads into a list to manipulate i.e. trim
        read_lst = []
        read_lst.append(read)

        new_lst = []

        for str in read_lst:
            for num_start in start_lst:
                for num_end in end_lst:
                    s = str.replace(str[(num_start-1):num_end],"")
                    new_lst.append(s)

I do get the trimmed reads, but what I'm having trouble figuring out is how to write these trimmed reads with their corresponding headers to a new FASTA file.
Also if there is any way that I could improve my code or employ better strategy to approach the problem, I would love to know.
Best wishes.
Edit: I tried running the solution by the super helpful user @Steve.
The command I ran:
python3 tailTrimming.py --fasta tailtrim.fasta --output testtrim.trim.fasta

The command line error:
Traceback (most recent call last):
  File "/home/x/Desktop/fastq/tailTrimming.py", line 61, in <module>
    main()
  File "/home/x/Desktop/fastq/tailTrimming.py", line 58, in main
    SeqIO.write(trim_reads(fasta), out, "fasta")
  File "/home/x/.local/lib/python3.10/site-packages/Bio/SeqIO/__init__.py", line 478, in write
    count = writer_class(fp).write_file(sequences)
  File "/home/x/.local/lib/python3.10/site-packages/Bio/SeqIO/Interfaces.py", line 209, in write_file
    count = self.write_records(records)
  File "/home/x/.local/lib/python3.10/site-packages/Bio/SeqIO/Interfaces.py", line 193, in write_records
    for record in records:
  File "/home/x/Desktop/fastq/tailTrimming.py", line 20, in trim_reads
    result = match.groupdict()
AttributeError: 'NoneType' object has no attribute 'groupdict'

This is the whole input file I supplied:
>SRR5655563.9676 9676 length=126 (type=T,start=1,end=11,length=11,identity=81.8182%)
tttctttctttAGTTCCAATAGTCTTGACCACCTTTTATCAACAGAAGCTTGGCCTTGTCCAGTTATATATTGTATTGTACCCTCAGATTTCTCGCACAAT
>SRR5655563.10087 10087 length=126 (type=T,start=2,end=16,length=15,identity=80%)
GtttgtttcattttttCCAAGGAAGAAATTCGACCTTCCAGCGATTTGTGGAGGTCATTCATCTGTTGAGTAAAGGATTGGTGCGCGACCTTGTAATCCTTGTGGAGGGTGTAG
>SRR5655563.13902 13902 length=126 (type=T,start=2,end=16,length=15,identity=80%)
CtttttcttcttgtttGATGTAAAAAATCTTCGGTATCTAAAACGTATGAAACCTTAGATGAACTGCAAGCTATCTTACACTCTAGACCATTGATAATTT
>SRR5655563.18489 18489 length=126 (type=T,start=2,end=12,length=11,identity=81.8182%)
CtttttcttcttCTTGGAAGTTTTGGATTCTGCAGGGACATTTTCTTTCTTTTCTCCATCTCCATCAGTCTCAACACCTGTCTTATCCTGAGCACCTTGAGACTCATCTTGACCATCTTTACTTT
>SRR5655563.20782 20782 length=126 (type=T,start=2,end=15,length=14,identity=71.4286%)
CttcattttagttttCAACTCCTCTACCAACTGGGATGGTGCATTCAACTGTAGAACTGTTTCACGAATCTCCTTGAGAACACTGAAGTCTCCCTCAGTTAACTTCTTTTTCAAATTATTGA
>SRR5655563.23972 23972 length=126 (type=T,start=2,end=29,length=28,identity=75%)
CttttttgttgttctcttctttcttctttGCATGGCTTGATGATGACTGATCCGATTGATTATTTCTTGATAGAATATCAATGTCATCACCAATTGTTGAGT
>SRR5655563.31739 31739 length=126 (type=T,start=1,end=10,length=10,identity=80%)
ttttcttcttCTTGCATGAAAGTGGGTTTGGACCCTTAGCTCTTTTTCTCTTGAACTGAGGTTT
>SRR5655563.34244 34244 length=126 (type=T,start=2,end=14,length=13,identity=76.9231%)
CtttcctttattttCTTGGAGGTTAGCCAGGATATTGGGATCAAAGCATGTAGCAAAACACACCTTGCTATTTCTTGCTCGTACAGATACATGGGATAGGA
>SRR5655563.38334 38334 length=126 (type=T,start=2,end=12,length=11,identity=81.8182%)
CtttattgttttCTCATCACCCATATTATATTTCTTAAGGGCAGATTCCAACTCTTCCATTGTGATGTACCCGCTCTTATCATTGTCAAAATATTCAAAGGCTTTATATAGGCGATCCTCTCTTTC

And this is the whole output I have gotten:
>SRR5655563.9676 9676 length=126 (type=T,start=1,end=11,length=11,identity=81.8182%)
AGTTCCAATAGTCTTGACCACCTTTTATCAACAGAAGCTTGGCCTTGTCCAGTTATATAT
TGTATTGTACCCTCAGATTTCTCGCACAAT

Just one read, that's all. Whereas I would like to keep all the sequences because all of them have this start, end attribute. The sequence being on two lines in the output is not something I did. It is like this in the output file. So it's no longer FASTA file, right?
Thanks in advance for any help!

I think there's a typo in your example. The part to trim does not exist in the input sequence. Also, if you want to keep the remaining sequence (i.e. "the rest of it"), shouldn't the output include the sequence before the starting position? Or would you like to extract the sequence after the end position?

Just to add in terms of performance Python is slow (loops and stuff). re is slow even for Python. Although Perl is not popular it will outperform the solutions here on a single CPU. If there are performance issues it's worth posting back - ultimately it depends on the size of the file and sensitivity to runtime. Huge files can simply fall over (i.e. never complete). However, I agree with Steve's solution: what Python is very cool at is splicing a string and that is super important to the solution. This is the bit that is really missing in the above code IMO

You made mistake describing the problem (sequence does not match 2-12 position intervals). If you update it, maybe you get more help. Also, it's not clear what kind of header you want to print in the record, exact same FASTA header?

i have the same problem but looks like i found some interesting solution @Steve. However, i have want to ask... what does this mean? pattern = ( r'type=(?P\w+),' r'start=(?P\d+),' r'end=(?P\d+),' r'length=(?P\d+),' r'identity=(?P\d+(?:.\d+)?)%' )
r'type= ??? what should i put here? r'start = ??? refers to the forward primer?
r'end = reverse primer?
r'length = the length of the final sequence after the trimming? what happens if there are sequences with different length afterwards? r'identity = the similarity between the primer and its sequence?

Here's one way using Biopython and the SeqIO interface to read and write SeqRecord objects. Using a generator, we can produce (yield) trimmed sequence records that can then be written to the output handle. This avoids storing an entire list of trimmed sequences in memory and is usually the way to go, assuming you're dealing with large files containing thousands of records:
import argparse
import re

from Bio import SeqIO

def trim_reads(fileobj):

    pattern = (
        r'type=(?P<type>\w+),'
        r'start=(?P<start>\d+),'
        r'end=(?P<end>\d+),'
        r'length=(?P<length>\d+),'
        r'identity=(?P<identity>\d+(?:\.\d+)?)%'
    )

    for seq_record in SeqIO.parse(fileobj, "fasta"):
        match = re.search(pattern, seq_record.description)
        result = match.groupdict()

        start = int(result['start']) - 1
        end = int(result['end'])

        head = seq_record.seq[:start]
        tail = seq_record.seq[end:]

        seq_record.seq = head + tail

        yield seq_record

def get_argument_parser():

    parser = argparse.ArgumentParser()

    parser.add_argument(
        "--fasta",
         type=str,
         required=True,
         help="input fasta file",
    )
    parser.add_argument(
        "--output",
        type=str,
        required=True,
        help="output tail.trimmed.fasta file",
    )

    return parser

def main():

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

    with open(args.fasta) as fasta, open(args.output, "w") as out:
        SeqIO.write(trim_reads(fasta), out, "fasta-2line")

if __name__ == "__main__":
    main()

Its worth point out that re is notoriously slow and over large files parallelisation is needed across the Biopython object.

Hello Steve. Thank you so much for you reply but I tried running your code and it gave the following error: AttributeError: 'NoneType' object has no attribute 'groupdict'

No worries at all, @pubsurfted. Yes - re.search will return None if the search pattern could not be found. You'll need to add some code to handle these cases. You might need a different search pattern, or maybe you just want to ignore these records? If, in your question above, you can provide a few examples of where the search pattern fails, I might be able to suggest a solution.

Hello @Steve, I have edited the question. If you could please take a look, thank you.

@pubsurfted I've updated the regex above to handle the fact that the identity can also be a plain integer, i.e. the decimal places are optional. Also, FASTA files are often line wrapped at 60 characters, but this entirely optional. These are sometimes called multiline FASTA files. There's probably a way to change the line wrapping using Biopython.

@pubsurfted There is indeed, using the "fasta-2line" format.

Thank you so so much Steve! Also, the "fasta-2line" format is raising the error: ValueError: Unknown format 'fasta-2line'

It ran. I had to update the biopython package from 1.67 to 1.80.