text:
Get frequency of index sequences of fastq file
I recently did some sequencing on a MiniSeq and unexpectedly the Undetermined_S0_L001_R1_001.fastq.gz file contained quite a lot of data (about 40% of total). I would like to know which indices are present in this file and ideally get a frequency of all index sequences in this file. How can I do this?
This is the output of head:
@MN00148:128:000H2YCW2:1:11102:12116:1065 1:N:0:0
GNAGAGGATTCAGGGAGCCAGCGCGCATATATCAG
+
F#FFFFFFFFFFFFFFFFFFA=/A///////////
@MN00148:128:000H2YCW2:1:11102:11354:1065 1:N:0:0
GNACCCGTGCGGCTGGAGAGACCACTAGAAAGAAT
+
A#AFFFFFFFFFFFFFFFFF///////FF=////6
@MN00148:128:000H2YCW2:1:11102:4103:1066 1:N:0:0
ANGTTGGTGGTACATCAGGGACGCGGAGATCTCAG

How did you demultiplex? With bcl2fastq? You can generate the index reads like I1

I actually went back to the output folder on the MiniSeq and found this file: DemultiplexSummaryF1L1, which contains demultiplexing information of the whole run. Half-way through this file I found what I was looking for:
### Most Popular Index Sequences

### Columns: Sequence ReverseComplement HitCount

Index
TGACCA  TGGTCA  2455964
CGATGT  ACATCG  2320579
GGCCCC  GGGGCC  2023652
CGCGGG  CCCGCG  1983605
CCGTCC  GGACGG  1641274
CAGATC  GATCTG  1514252
CCGGCC  GGCCGG  1319383
CCGCGC  GCGCGG  1109940
ACAGTG  CACTGT  891344
CCCGGG  CCCGGG  706083
TGCCCA  TGGGCA  168294
CAGCTC  GAGCTG  164091
.... [data continues]

I don't know that you can get them from that file.  I've noticed that when the Miseq is left to run its version of fastq generation software, there are no indices in the undetermined file.  I usually run bcl2fastq myself if I want that information.
They recently came out with a new version of software for the MiSeq, so things might be different now.


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

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

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

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

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

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

source: Stack Exchange

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

sort: votes


text:
Estimate insert size for single-end data with picard CollectInsertSizeMetrics
I have a BAM file generated from single-end data and I want to estimate the insert size using picard's CollectInsertSizeMetrics as follows:
java -jar Picard/picard.jar CollectInsertSizeMetrics I=Ctrl_sorted.bam O=Ctrl.insert_sizes.txt H=Ctrl.insert_sizes.pdf M=0.5 MAX_RECORDS_IN_RAM=2000000 VALIDATION_STRINGENCY=LENIENT

It runs fine until at the very end I get the following message:
WARNING 2023-01-31 15:37:13 CollectInsertSizeMetrics    All data categories were discarded because they contained < 0.5 of the total aligned paired data.

Picard obviously thinks the data should be paired-end, so is there a way to tell it that it is single-end? I cannot find an option for this:
java -jar picard.jar CollectInsertSizeMetrics -h
USAGE: CollectInsertSizeMetrics [arguments]

Collect metrics about the insert size distribution of a paired-end library. This tool provides useful metrics for
validating library construction including the insert size distribution and read orientation of paired-end
libraries.
The expected proportions of these metrics vary depending on the type of library preparation used,
resulting from technical differences between pair-end libraries and mate-pair libraries. For a brief primer on
paired-end sequencing and mate-pair reads, see the GATK Dictionary.
The CollectInsertSizeMetrics
tool outputs the percentages of read pairs in each of the three orientations (FR, RF, and TANDEM) as a histogram. In
addition, the insert size distribution is output as both a histogram (.insert_size_Histogram.pdf) and as a data table
(.insert_size_metrics.txt).

Note: Metrics labeled as percentages are actually expressed as fractions!

Usage
example:

java -jar picard.jar CollectInsertSizeMetrics \
      I=input.bam \

O=insert_size_metrics.txt \
      H=insert_size_histogram.pdf \
      M=0.5
Note: If processing a small
file, set the minimum percentage option (M) to 0.5, otherwise an error may occur.

Please see InsertSizeMetrics
for detailed explanations of each metric.

Version:2.23.8

Required Arguments:

--Histogram_FILE,-H:File      File to write insert size Histogram chart to.  Required.

--INPUT,-I:File               Input SAM or BAM file.  Required.

--OUTPUT,-O:File              The file to write the output to.  Required.

Optional Arguments:

--arguments_file:File         read one or more arguments files and add them to the command line  This argument may be
                              specified 0 or more times. Default value: null.

--ASSUME_SORTED,-AS:Boolean   If true (default), then the sort order in the header file will be ignored.  Default value:
                              true. Possible values: {true, false}

--COMPRESSION_LEVEL:Integer   Compression level for all compressed files created (e.g. BAM and VCF).  Default value: 5.

--CREATE_INDEX:Boolean        Whether to create an index when writing VCF or coordinate sorted BAM output.  Default
                              value: false. Possible values: {true, false}

--CREATE_MD5_FILE:Boolean     Whether to create an MD5 digest for any BAM or FASTQ files created.    Default value:
                              false. Possible values: {true, false}

--DEVIATIONS:Double           Generate mean, sd and plots by trimming the data down to MEDIAN +
                              DEVIATIONS*MEDIAN_ABSOLUTE_DEVIATION. This is done because insert size data typically
                              includes enough anomalous values from chimeras and other artifacts to make the mean and sd
                              grossly misleading regarding the real distribution.  Default value: 10.0.

--GA4GH_CLIENT_SECRETS:String Google Genomics API client_secrets.json file path.  Default value: client_secrets.json.

--help,-h:Boolean             display the help message  Default value: false. Possible values: {true, false}

--HISTOGRAM_WIDTH,-W:Integer  Explicitly sets the Histogram width, overriding automatic truncation of Histogram tail.
                              Also, when calculating mean and standard deviation, only bins <= Histogram_WIDTH will be
                              included.  Default value: null.

--INCLUDE_DUPLICATES:Boolean  If true, also include reads marked as duplicates in the insert size histogram.  Default
                              value: false. Possible values: {true, false}

--MAX_RECORDS_IN_RAM:Integer  When writing files that need to be sorted, this will specify the number of records stored
                              in RAM before spilling to disk. Increasing this number reduces the number of file handles
                              needed to sort the file, and increases the amount of RAM needed.  Default value: 500000.

--METRIC_ACCUMULATION_LEVEL,-LEVEL:MetricAccumulationLevel
                              The level(s) at which to accumulate metrics.    This argument may be specified 0 or more
                              times. Default value: [ALL_READS]. Possible values: {ALL_READS, SAMPLE, LIBRARY,
                              READ_GROUP}

--MIN_HISTOGRAM_WIDTH,-MW:Integer
                              Minimum width of histogram plots. In the case when the histogram would otherwise
                              betruncated to a shorter range of sizes, the MIN_HISTOGRAM_WIDTH will enforce a minimum
                              range.  Default value: null.

--MINIMUM_PCT,-M:Float        When generating the Histogram, discard any data categories (out of FR, TANDEM, RF) that
                              have fewer than this percentage of overall reads. (Range: 0 to 1).  Default value: 0.05.

--QUIET:Boolean               Whether to suppress job-summary info on System.err.  Default value: false. Possible
                              values: {true, false}

--REFERENCE_SEQUENCE,-R:File  Reference sequence file.  Default value: null.

--STOP_AFTER:Long             Stop after processing N reads, mainly for debugging.  Default value: 0.

--TMP_DIR:File                One or more directories with space available to be used by this program for temporary
                              storage of working files  This argument may be specified 0 or more times. Default value:
                              null.

--USE_JDK_DEFLATER,-use_jdk_deflater:Boolean
                              Use the JDK Deflater instead of the Intel Deflater for writing compressed output  Default
                              value: false. Possible values: {true, false}

--USE_JDK_INFLATER,-use_jdk_inflater:Boolean
                              Use the JDK Inflater instead of the Intel Inflater for reading compressed input  Default
                              value: false. Possible values: {true, false}

--VALIDATION_STRINGENCY:ValidationStringency
                              Validation stringency for all SAM files read by this program.  Setting stringency to
                              SILENT can improve performance when processing a BAM file in which variable-length data
                              (read, qualities, tags) do not otherwise need to be decoded.  Default value: STRICT.
                              Possible values: {STRICT, LENIENT, SILENT}

--VERBOSITY:LogLevel          Control verbosity of logging.  Default value: INFO. Possible values: {ERROR, WARNING,
                              INFO, DEBUG}

--version:Boolean             display the version number for this tool  Default value: false. Possible values: {true,
                              false}

Advanced Arguments:

--showHidden,-showHidden:Boolean
                              display hidden arguments  Default value: false. Possible values: {true, false}

Or is there another tool that can deal with single-end data to estimate the insert size?

I can't see how you would calculate insert size from single end data. Are you aware of a paper about it, or something like that?

I just realised I wasn't thinking clearly; you are right - from single-end data you cannot calculate the insert size. This can only be established experimentally. Thanks everyone for your comments.

CollectInsertSizeMetrics does not estimate but (as by the name) collects insert size metrics, which is nothing different than parsing the TLEN field from the BAM file. It is intended to only be used for paired-end data, not single-end, as single-end data do not have a value in that field in the BAM file. It is calculated by the aligner for paired-end data.
What kind of data do you have? For targeted assays such as ChIP-seq the macs2 peak caller has a module to estimate fragment size based on strand cross correlation analysis. I do not kow whether you can use that for any sort of assay.
If this is your own data then you can ask the lab person who made the library to look up the size from the Bioanalyzer or TapeStation QC. If it is published data then check the methods text for the procedure/kit. Most NGS types have a relatively standard insert size depending on the exact assay.

It is from published data (RNA-Seq). I want to do alternative splicing analysis with MISO, which asks for the insert size.
Thanks for your suggestions, I will look deeper into the methods.

Can you link the GEO reference? RNA-seq libraries usually have a somewhatish 200-300bp insert size.

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE121949

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


text:
Most reads no_feature with htseq-count
I have mapped reads to the yeast (R64-1-1) genome with STAR like this:
STAR --runThreadN 44 --runMode alignReads --genomeDir /home/user/Documents/scripts/pyseqtools/index/STAR/R64-1-1 --readFilesIn /home/user/Documents/analyses/TT-Seq/user/trim/WT_3_R1_001_val_1.fq.gz /home/user/Documents/analyses/TT-Seq/user/trim/WT_3_R2_001_val_2.fq.gz --readFilesCommand zcat --quantMode TranscriptomeSAM GeneCounts --twopassMode Basic --outSAMunmapped None --outSAMattrRGline ID:WT_3 PU:WT_3 SM:WT_3 LB:unknown PL:illumina --outSAMtype BAM Unsorted --outTmpDir /home/user/Documents/analyses/TT-Seq/user/tmp --outFileNamePrefix /home/user/Documents/analyses/TT-Seq/user/bam/R64-1-1/WT_3/WT_3

I wanted to count the reads using htseq-count as follows:
htseq-count --quiet --format bam --stranded yes /home/user/Documents/analyses/TT-Seq/user/bam/R64-1-1/WT_3/WT_3Aligned.out.bam /home/user/Documents/scripts/pyseqtools/gtf/R64-1-1/Saccharomyces_cerevisiae.R64-1-1.106.gtf > /home/user/Documents/analyses/TT-Seq/user/htseq-count/WT_3.count.txt

The problem is that most reads map to no_feature:
__no_feature    4462786
__ambiguous 12644
__too_low_aQual 0
__not_aligned   4907
__alignment_not_unique  12529444

I have check whether the chromosome names in the BAM file and the GTF file matched (which they seem to do):
BAM:
@HD VN:1.4
@SQ SN:I    LN:230218
@SQ SN:II   LN:813184
@SQ SN:III  LN:316620
@SQ SN:IV   LN:1531933
@SQ SN:V    LN:576874
@SQ SN:VI   LN:270161
@SQ SN:VII  LN:1090940
@SQ SN:VIII LN:562643
@SQ SN:IX   LN:439888
@SQ SN:X    LN:745751
@SQ SN:XI   LN:666816
@SQ SN:XII  LN:1078177
@SQ SN:XIII LN:924431
@SQ SN:XIV  LN:784333
@SQ SN:XV   LN:1091291
@SQ SN:XVI  LN:948066
@SQ SN:Mito LN:85779

GTF:
IV  sgd gene    8683    9756    .   -   .   gene_id "YDL246C"; gene_name "SOR2"; gene_source "sgd"; gene_biotype "protein_coding";
IV  sgd transcript  8683    9756    .   -   .   gene_id "YDL246C"; transcript_id "YDL246C_mRNA"; gene_name "SOR2"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "SOR2"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd exon    8683    9756    .   -   .   gene_id "YDL246C"; transcript_id "YDL246C_mRNA"; exon_number "1"; gene_name "SOR2"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "SOR2"; transcript_source "sgd"; transcript_biotype "protein_coding"; exon_id "YDL246C_mRNA-E1";
IV  sgd CDS 8686    9756    .   -   0   gene_id "YDL246C"; transcript_id "YDL246C_mRNA"; exon_number "1"; gene_name "SOR2"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "SOR2"; transcript_source "sgd"; transcript_biotype "protein_coding"; protein_id "YDL246C";
IV  sgd start_codon 9754    9756    .   -   0   gene_id "YDL246C"; transcript_id "YDL246C_mRNA"; exon_number "1"; gene_name "SOR2"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "SOR2"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd stop_codon  8683    8685    .   -   0   gene_id "YDL246C"; transcript_id "YDL246C_mRNA"; exon_number "1"; gene_name "SOR2"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "SOR2"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd gene    17577   18566   .   -   .   gene_id "YDL243C"; gene_name "AAD4"; gene_source "sgd"; gene_biotype "protein_coding";
IV  sgd transcript  17577   18566   .   -   .   gene_id "YDL243C"; transcript_id "YDL243C_mRNA"; gene_name "AAD4"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "AAD4"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd exon    17577   18566   .   -   .   gene_id "YDL243C"; transcript_id "YDL243C_mRNA"; exon_number "1"; gene_name "AAD4"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "AAD4"; transcript_source "sgd"; transcript_biotype "protein_coding"; exon_id "YDL243C_mRNA-E1";
IV  sgd CDS 17580   18566   .   -   0   gene_id "YDL243C"; transcript_id "YDL243C_mRNA"; exon_number "1"; gene_name "AAD4"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "AAD4"; transcript_source "sgd"; transcript_biotype "protein_coding"; protein_id "YDL243C";
IV  sgd start_codon 18564   18566   .   -   0   gene_id "YDL243C"; transcript_id "YDL243C_mRNA"; exon_number "1"; gene_name "AAD4"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "AAD4"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd stop_codon  17577   17579   .   -   0   gene_id "YDL243C"; transcript_id "YDL243C_mRNA"; exon_number "1"; gene_name "AAD4"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "AAD4"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd gene    1248154 1249821 .   -   .   gene_id "YDR387C"; gene_name "CIN10"; gene_source "sgd"; gene_biotype "protein_coding";
IV  sgd transcript  1248154 1249821 .   -   .   gene_id "YDR387C"; transcript_id "YDR387C_mRNA"; gene_name "CIN10"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "CIN10"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd exon    1248154 1249821 .   -   .   gene_id "YDR387C"; transcript_id "YDR387C_mRNA"; exon_number "1"; gene_name "CIN10"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "CIN10"; transcript_source "sgd"; transcript_biotype "protein_coding"; exon_id "YDR387C_mRNA-E1";
IV  sgd CDS 1248157 1249821 .   -   0   gene_id "YDR387C"; transcript_id "YDR387C_mRNA"; exon_number "1"; gene_name "CIN10"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "CIN10"; transcript_source "sgd"; transcript_biotype "protein_coding"; protein_id "YDR387C";
IV  sgd start_codon 1249819 1249821 .   -   0   gene_id "YDR387C"; transcript_id "YDR387C_mRNA"; exon_number "1"; gene_name "CIN10"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "CIN10"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd stop_codon  1248154 1248156 .   -   0   gene_id "YDR387C"; transcript_id "YDR387C_mRNA"; exon_number "1"; gene_name "CIN10"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "CIN10"; transcript_source "sgd"; transcript_biotype "protein_coding";

Edit: the output of ReadsPerGene.out.tab is (thanks
swbarnes2) :
N_unmapped  73359652    73359652    73359652
N_multimapping  15789900    15789900    15789900
N_noFeature 192619  5872967 380420
N_ambiguous 648956  13479   430527

Would anyone know why most reads map to no_feature?

You are already running it with GeneCounts, which should output the same as htseq-count.  What does your ReadsPerGene.out.tab say?

Do your gtf contains exon entris (i.e. with "exon" on the 3rd column)? Could you please update your question with a few exon entries from your gtf?


text:
Subset GTF file for specific genes
I have a GTF file (in 'test.gtf' file):
chr1    HAVANA  gene    10069   18412   .   +   .   gene_id "ENSG00000223972.4"; transcript_id "ENSG00000223972.4"; gene_type "pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1"; level 2; havana_gene "OTTHUMG00000000961.2";
chr2    HAVANA  gene    11869   14412   .   +   .   gene_id "ENSG0000023172.1"; transcript_id "ENSG000023172.1"; gene_type "pseudogene"; gene_status "KNOWN"; gene_name "PCNA"; transcript_type "pseudogene"; transcript_status "KNOWN"; transcript_name "PCNA"; level 2; havana_gene "OTTHUMG00000001931.2";
chr3    HAVANA  gene    15769   19412   .   +   .   gene_id "ENSG00000433972.1"; transcript_id "ENSG000023172.1"; gene_type "pseudogene"; gene_status "KNOWN"; gene_name "USP21"; transcript_type "pseudogene"; transcript_status "KNOWN"; transcript_name "USP21"; level 2; havana_gene "OTTHUMG00000002361.2";

I have loaded it as a pyranges object and I want to subset it so that it only contains gene from test_list:
import pyranges as pr
import pandas as pd

gtf_file = "test.gtf"
test_list = ["PCNA","USP21", "USP1"]

gtf = pr.read_gtf(gtf_file)

sub_gtf = gtf

sub_gtf.df = gtf.df[gtf.df.gene_name.isin(test_list)]

But the last line throws this error:
File "<ipython-input-24-af0b4b5c7576>", line 1, in <module>
    irongtf.df = irongtf.df[irongtf.df.gene_name.isin(test_list)]

  File "/home/niek/.local/lib/python3.8/site-packages/pyranges/pyranges.py", line 314, in __setattr__
    _setattr(self, column_name, column)

  File "/home/niek/.local/lib/python3.8/site-packages/pyranges/methods/attr.py", line 48, in _setattr
    df.insert(pos, column_name, column)

  File "/home/niek/.local/lib/python3.8/site-packages/pandas/core/frame.py", line 4413, in insert
    value = self._sanitize_column(value)

  File "/home/niek/.local/lib/python3.8/site-packages/pandas/core/frame.py", line 4504, in _sanitize_column
    com.require_length_match(value, self.index)

  File "/home/niek/.local/lib/python3.8/site-packages/pandas/core/common.py", line 531, in require_length_match
    raise ValueError(

ValueError: Length of values (2) does not match length of index (1)

What would be the right way to subset the gtf file so that it only contains the rows where gene_name contains the gene symbols from test_list?
Or is there another tool that I could use to do this?

It turned out to be a lot easier if I used grep to do this:
grep -w -f genes.txt gencode.v19.annotation.gtf > sub_set.gtf

genes.txt contains gene symbols on new lines.

the only risk you are taking with this approach is that some of the gene names might be sub-strings of other gene names and as a consequence, those would get pulled as well.

I thought that the -w flag would prevent that?

Ah, you are right, I overlooked -w. Should be fine, unless the genes names contain non-word constituent characters (like - or .).

Maybe a better way would be: grep -wf <(awk '{print "gene_name \"$0\""}' genes.txt) gencode.v19.annotation,gtf > subset.gtf. This will ensure that the strings are compared only to the gene_name tag in each feature.

Great suggestion Ram RS; I have tried it but it gives me an empty file?

I'm not sure what the problem could be. Try just the awk command on the first 1000 lines of the big GTF, then try the grep with one of the output lines from the awk command. Step by step debugging like this is the way to solve these problems.

@justinian482 have you ever edited the genes.txt file in Windows? Could there be Windows line endings?

@terdon no I am using Ubuntu 20.04.3

Following up on zorbax's answer, you could read in and filter the GTF file in this way, among others:
#!/usr/bin/env python

import gtfparse as gp

gtf_file = "test.gtf"
test_list = ["PCNA", "USP21", "USP1"]

df = gp.read_gtf(gtf_file)
subset = df[df['gene_name'].str.contains('|'.join(test_list))]
print(subset)

You could further cull this subset by other attributes or columns, e.g.:
genes_only = subset[subset["feature"] == "gene"]

etc.
The gtfparse library is available via PyPi: https://pypi.org/project/gtfparse/ or Bioconda: https://bioconda.github.io/recipes/gtfparse/README.html

How would I save the df as a proper GTF after subsetting? The separators of the final columns are different compared to the first ones (; vs \t)?

You could add import pandas as pd and then try df.to_csv(out_filename, sep='\t') to write out a tab-delimited file from various data frame columns. You'll probably need to first reduce the attribute columns to a single column (https://stackoverflow.com/questions/27298178/concatenate-strings-from-several-rows-using-pandas-groupby), and then remove those additional columns. Or create a new data frame from the non-attribute and attribute columns and export that. Perhaps try it out and if you get stuck, ask another question showing what you've tried.

If you have your gtf file like a DataFrame you can use:
 df[df['gene_name'].str.contains('|'.join(test_list))]


text:
Filtering samples by FORMAT attributes using bcftools
The -e and -i options of the bcftools filter command appear, by default, to only allow for including or excluding sites.  E.g., -e 'FMT/DP < 10' removes sites where any sample has DP < 10, and -e 'MEAN(FMT/DP) < 10' removes sites where average depth across samples is < 10.
I would like to perform effectively similar filtering commands, but in a way that includes or excludes samples, instead of sites.  E.g., exclude all samples that have DP < 10 at any site, or exclude all samples with an average depth across all sites < 10, or exclude all samples that have a genotype quality of < 20 at 10% of genotyped sites.
Such commands can very easily be written to apply to sites, but there does not seem to be a straightforward way of applying such filtering expressions to samples in bcftools.
One option mentioned in the manual is to use the double || operator, as described in the man page.  As described in my comment on a Github issue, however, this does not seem to work as a filtering operation for including or excluding entire columns in a VCF file, but rather only for selecting more or less individual records (e.g., as illustrated in the example on this tutorial page).  In addition, it would seem limited in functionality (e.g., functions such as MEAN(FMT/DP) < 10 would still be evaluated across samples, not sites, when combined with a || operator).
The only other method mentioned in the manual is to use the SMPL or s prefix to functions which evaluate FORMAT fields in expressions passed to filter.  This also does not appear to work as I would expect.  I'm not sure I understand what it is doing, but it definitely isn't removing columns from the VCF file.  Instead it seems to be using as a metric calculated across sites on a per sample basis as a filtering criteria applied to sites (but I could just be confused):
nSites_unfiltered=$(bcftools view -H unfiltered.vcf.gz | wc -l)
nSamples_unfiltered=$(bcftools query unfiltered.bcf.gz -l | wc -l)

nSites_DPfilter=$(bcftools filter unfiltered.bcf.gz -e 'sMEAN(FORMAT/DP)<5' | bcftools view -H | wc -l)
nSamples_DPfilter=$(bcftools filter unfiltered.bcf.gz -e 'sMEAN(FORMAT/DP)<5' | bcftools query -l | wc -l)

echo "scale=5; $nSites_DPfilter / $nSites_unfiltered" | bc #0.17349 - sites are removed, even though the filter is being applied to samples
echo "scale=5; $nSamples_DPfilter / $nSamples_unfiltered" | bc #1.00000 - no samples removed, despite, again, the filter being applied to samples

nSites_normalDPFilter=$(bcftools filter unfiltered.bcf.gz -e 'MEAN(FMT/DP) < 5' | bcftools view -H | wc -l)
echo "scale=5; $nSites_normalDPFilter / $nSites_unfiltered" | bc #0.83121 - fewer sites removed than the sMEAN filter

The .bcf.gz passed to these commands is here: https://www.dropbox.com/s/fwj1snabsb763sn/unfiltered.bcf.gz?dl=0
Its .csi index file is here: https://www.dropbox.com/s/jzdq6y91nkencmy/unfiltered.bcf.gz.csi?dl=0
In any event, the only option I have been able to come up with that actually excludes or includes entire sample IDs is to extract summary statistics using bcftools stats, export the desired statistics from the PSC, Per-sample counts table, and perform my own exploratory analysis outside of bcftools.  Then I can use bcftools view -S to explicitly remove samples by their ID, which I have identified manually as failing some thresholding criteria.  This also is obviously cumbersome, and only works for FMT fields that are included in stats (e.g., not GQ or PL).
I am not sure if I am missing some basic functionality of bcftools, or if this is more of a feature request.  Any information would be appreciated.
Thanks!

You may want to take a look at vcftools' genotype filtering options or the custom JEXL expressions that GATK's SelectVariants allows. Are you sure your sample-filtering operations won't produce a variable number of columns for each row? Also, remember that bcftools filter will not remove samples no matter what you try. It's a utility to filter sites, and it can only filter sites.

Thanks very much.  Since every site has genotypes for every sample (if you include missing genotype entries), I'm not sure I understand how removing certain columns in their entirety would create any problems in terms of uneven dimensions.  Indeed, this is the function of bcftools view -s ^<list.of.sample.IDs>, if I am not mistaken.  vcftools also only allows for filtering samples based on ID.  Thank you for the reference to SelectVariants.  Do you have a link to a manual, or example code to address any of the example use-cases referenced in my question?  Thanks!

Not really, no. I have rarely ever used SelectVariants when it wasn't already implemented in a pipeline (GATK tools can be a pain to use at times). If push comes to shove, it might be worth converting the VCF to MAF format so the data is at individual level.


text:
Quantify low complexity regions in DNA sequences
I have fasta files with multiple sequences. They are reads that mapped to the genome outside of  probe-targeted regions.
From a quick perusal, they appear to be repetitive and have low complexity. Is there a way I can quantify their complexity and use that as an explanation as to why those regions were mapped to?

Use Picard EstimateLibraryComplexity.
From the docs:

Estimates the numbers of unique molecules in a sequencing library.

java -jar picard.jar EstimateLibraryComplexity \
     I=input.bam \
     O=est_lib_complex_metrics.txt

Although I mentioned I had fasta files, I was able to take the original bam files and pipe only those regions into picard with ´samtools view -b -L regions.bed sample.bed | picard....´

"Complexity" is also known as entropy. A quick and dirty way is to use a file compression utility - zip will get good compression for a repetitive sequence.
You can also look at the k-mer distribution for some suitable k: in a repetitive sequence it will be biased. Here is a tool that give a local plot: https://rdrr.io/github/vsbuffalo/qrqc/man/kmerEntropyPlot.html

I looked into compressibility as a measure of complexity. I found this article very interesting. They do not recommend it for DNA sequences.
https://sci-hub.do/10.1109/iwobi.2014.6913941

How long are your reads, and how many?  Some colleagues of mine wrote https://github.com/eclarke/komplexity to quickly score complexity in large numbers of short reads.  The score is derived from the number of unique k-mers normalized by the sequence length, so it's simple and fast, though the score inherently becomes less informative for longer and longer sequences (so, best for short reads).  We've used it for filtering repetitive or other low-complexity sequences in metagenomics projects.

This looks promising!
I don't want to run it on the unmapped per se as you can see in my question. But since all of the non-targeted regions are less than 300 bps, I think komplixity will be helpful.  I'll let  you know how it goes, and possibly award this as the accepted answer.

do you think you could ask your colleagues to take a look at their issues on github? I'm not sure it's being maintained.

I see what you mean, those issues aren't getting any attention.  I think the short answer for most of those is that it'd need new features to address them, like handling paired-end reads automatically.  (But for that use case you can at least filter forward and reverse separately, and take the intersection of what remains.)  If I can carve out the time I'll send eclarke some PRs to help with those though.

Thanks Jesse. That would be helpful.


text:
Use biopython to align SeqRecords stored in dict
I'd like to perform multiple alignments, where a gene from each sample was read in from fasta files. The fasta file represented one sample and had multiple genes. I have read in each sample fasta file and now have a dictionary of genes and their samples and sequences.
Here is a small example.
The main key is the gene ID for another dictionary as the value. In the dictionary 'R??' is the sample name and key,
followed by the annotation and sequence in a list.
{'PRELSG_MIT_v1:5894-6032': {'R47': ['not annotated',
                                     Seq('ATAAAAATAGTTACCATAGCTGTTGATGGATGCTTCGATCATAAAATATGTTTG...GGA')],
                             'R48': ['not annotated',
                                     Seq('ATAAAAATAGTTACCATAGCTGTTGATGGATGCTTCGATCATAAAATATGTTTG...GGA')],
                             'R49': ['not annotated',
                                     Seq('ATAAAAATAGTTACCATAGCTGTTGATGGATGCTTCGATCATAAAATATGTTTG...GGA')],
                             'R50': ['not annotated',
                                     Seq('ATAAAAATAGTTACCATAGCTGTTGATGGATGCTTCGATCATAAAATATGTTTG...GGA')],
                             'R51': ['not annotated',
                                     Seq('ATAAAAATAGTTACCATAGCTGTTGATGGATGCTTCGATCATAAAATATGTTTG...GGA')],
                             'R52': ['not annotated',
                                     Seq('ATAAAAATAGTTACCATAGCTGTTGATGGATGCTTCGATCATAAAATATGTTTG...GGA')]},
 'exon_PRELSG_MIT00300.1-E1': {'R47': ['cytochrome b',
                                       Seq('TCCAGATAATGCTATTACAGTAGATAGATATGCTACTCCTTTACATATTGTTCC...TTC')],
                               'R48': ['cytochrome b',
                                       Seq('TCCAGATAATGCTATTACAGTAGATAGATATGCTACTCCTTTACATATTGTTCC...TTC')],
                               'R49': ['cytochrome b',
                                       Seq('TCCAGATAATGCTATTACAGTAGATAGATATGCTACTCCTTTACATATTGTTCC...TTC')],
                               'R50': ['cytochrome b',
                                       Seq('TCCAGATAATGCTATTACAGTAGATAGATATGCTACTCCTTTACATATTGTTCC...TTC')],
                               'R51': ['cytochrome b',
                                       Seq('TCCAGATAATGCTATTACAGTAGATAGATATGCTACTCCTTTACATATTGTTCC...TTC')],
                               'R52': ['cytochrome b',
                                       Seq('TCCAGATAATGCTATTACAGTAGATAGATATGCTACTCCTTTACATATTGTTCC...TTC')]}}

I'd like to align each gene and write out those alignments to file without any intermediate files written to disk.
How do I do that?

do you want to align each of those Seq('ATAAAAATAGTT...... on top of each other ?

Yes, what you said moth. I want one alignment for each gene.

I don't believe this is possible. So I'm answering this to close it, since I don't have enough reputation to close it myself.

This BioPython document should get you started.
A specific example is given here.

All the examples in documentation you linked require the files to be ready to be aligned, where each file represents only one gene and multiple samples. That is not what I have here, and I don't want to write intermediate gene files. I want to do it all in memory until the final alignments are ready.

Strange requirement. No idea how to align in python without calling another application by file.

I learnt how to create a MultipleSeqAlignment object based on that dict, but I can't see how to pass that object to a muscle subprocess...
I guess I'd need a muscle python library.

I don't think muscle can work 'file free'... may be wrong.

Yeah, I've given up on this. Biopython supports this using pairwise alignments, but that's obviously not going to work here. Thanks for your help.

Good, what you're trying is very weird.


text:
How to get a FASTA format file with the DNA sequences of all annotated genes?
I am analysing Pyrococcus Furiosus DNA sequencing data by considering data published here in NCBI. When I click on "Send to">"Gene Features">"FASTA format" I download a file that has the sequences of genes of this organism but I realized that this file has some sequences that are doubled (that is the sequence is present twice with different identifiers)... is there a way to get a file with all genes annotated and their respective DNA sequences without double sequences and with all sequences in the 5' to 3' direction without reverse complement in a well "ordered" way ? In NCBI (in the link I reported before) it indicates 2,128 genes so I would like a file with all these 2,128 genes annotated and their respective DNA sequences in FASTA format. Do you know if there is an other website or an other place in NCBI in which I can search to get this kind of file?
EDIT: if you follow the path "Send to">"Gene Features">"FASTA format", now it seems that there are just 2,128 genes so this problem is vanished in this case. Anyway, I would appreciate the answer for the future and for the reverse complement gene sequences that are still there.
Thank you very much.

From this SO post. You can use seqkit to remove duplicate sequences with the command below:
seqkit rmdup -s < sequences.txt > out.fa

The rmdup option removes duplicates, and the -s option calls duplicates on the basis of sequence, ignoring differences in headers.

thank you @Throckmorton I will try ! Even if my question was more focus on if it exists some website (in general) that has a FASTA format file with all annotated genes, without double DNA sequences or sequences that are reverse complement etc... so that if I remove the FASTA headers, I would obtain the complete genome DNA sequence

@Manuela does this species really have nothing but genes in its genome? And all of its genes are unique? In other words, are you sure that concatenating all genes will give you the genome? That seems very unlikely to me, but then I've spent most of my career working on eukaryotes.

@terdon thank you for the observation ! I assumed it was so... I am not a biologist :)

Ah, then @Manuela it is almost certainly the case the the genome has a lot of DNA that isn't part of any gene. I don't know archea at all, but if we take human as an example, genes are only around ~5% of the genome. If you want the genome, then look for the genome and not the genes.

ok @terdon thank you .

Sorry @terdon , but I have a doubt related to this: is it completely wrong to search a DNA sequence motif in a genome sequence  (because I could find a DNA sequence similar to my DNA sequence motif within the part that does not belong to gene and then I count it as a motif but it can not be) ?

@Manuela that depends entirely on what you are looking for, what species you are working on etc. If your motif has a function that is only relevant within a coding sequence, for example if you are looking for a DNA sequence which when translated produces a protein domain, then you want to focus on coding sequences only. If you want to see how common the motif is in general, then you can look at the whole genome. But there is no simple answer, it will always depend on what question you are asking.

ok thank you very much.


text:
How to modify dot plot in MUMmer 3 for bacteria comparative genomics?
I was trying to make a dotplot to visualize genome-genome sequence alignment by using the MUMmer 3 software (it is typically used to compare whole genome sequences of bacteria). I runned the same command lines that I wrote in this question post What is the file .mums created by mummer?.
I obtain a plot that is consistent with what I expected but I do not understand why it does not show the axis scales (or labels ? I do not know how to explain but I report below my plot vs the plot made in the MUMmer page example 2.1. mummer so that you can understand what I mean). There could be two solutions:
1. Since my command lines are the same as the one of the example that I found on web, I think that I should modify the mummerplot code to fix this problem.
However, the code is written in C/C++, so not really easy for me. I tried to modify some lines of code but with no success. If you want to have a look at the code and give me some suggestions, it is available on Github.
2. But I wonder also if it is the only way to adjust the mummerplot plot. I found in the MUMmer 3 manual this:

The mummerplot script outputs three files, .gp .fplot .rplot, when run with standard parameters. The first of which is the gnuplot script. This script contains the commands necessary to generate the plot, and refers to the two data files which contain the forward and reverse matches respectively.  [...] Line thickness, color, and many other options can be added or removed from the plot by hand editing the gnuplot script.

Indeed I noticed that my .gp file has some differences with respect to the .gp file of the example that I linked above (the one here). So I could simply modify my .gp file writing it as the one of the example and obtain the plot as I want it, as also the manual suggests. But how can I do it if running mummerplot to generate my plot I obtain the .ps image of the plot and the .gp file at the same time ? The .gp file is not something that I can create before running the mummerplot command...
Below I report, respectively, my plot and the plot of the example (whose aesthetics I would like to imitate):

Thus, I ask if you can give me any hint by answering to one of the two questions or even both or even if there is a third way to solve that I ignore. Thank you in advance !


text:
Help to create a dataframe in Python from a FASTA file
I want to create a dataframe in Python starting from a FASTA format file. Given the toy FASTA file that I am attaching, I built this program in Python that returns four colums corresponding to id, sequence length, sequence, animal name and rows corresponding to all the data available.
However, I am trying to understand how to modify this code in order to create a dataframe in which classes Human and Dog have the same number of data. For example, I want to say to Python: "Append to record (that is the empty list) id, sequence length, sequence and animal for Human, but do it a number of times that is equal to the number of data of the class with minimum number of data (that is Dog)".
I think that a while loop is needed but I have a bit troubles to understand how to do it. Any suggestion ?
Below the Python code I wrote and the FASTA format file I used.
import pandas as pd
import re
def read_fasta(file_path, columns) :
    from Bio.SeqIO.FastaIO import SimpleFastaParser
    with open("Proof.txt") as fasta_file :
        records = [] # create empty list
        for title, sequence in SimpleFastaParser(fasta_file): #SimpleFastaParser Iterate over Fasta records as string tuples. For each record a tuple of two strings is returned, the FASTA title line (without the leading ‘>’ character), and the sequence (with any whitespace removed).
            record = []
            title_splits=re.findall(r"[\w']+", title) # Data cleaning is needed

            record.append(title_splits[0])  #First values are ID (Append adds element to a list)
            record.append(len(sequence)) #Second values are sequences lengths
            sequence = " ".join(sequence) #It converts into one line
            record.append(sequence)#Third values are sequences

            #Fourth column will contain the species
            if "Human" in title_splits:
                    record.append("Human")
            else:
                    record.append("Dog")

            records.append(record)
    return pd.DataFrame(records, columns = columns) #We have created a function that returns a dataframe

#Now let's use this function by inserting in the first argument the file name (or file path if your working directory is different from where the fasta file is)
#And in the second one the names of columns
data = read_fasta("Proof.txt", columns=["id","sequence_length", "sequence", "animal"])
data

The FASTA format file is this:
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Dog|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Dog|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Dog|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Dog|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK

My code prints a dataframe like:
       id  sequence_length                               sequence animal
0   Numer               19  H S S F I E I V N I E H V I E H I V K  Human
1   Numer               19  H S S F I E I V N I E H V I E H I V K  Human
2   Numer               19  H S S F I E I V N I E H V I E H I V K  Human
3   Numer               19  H S S F I E I V N I E H V I E H I V K  Human
4   Numer               19  H S S F I E I V N I E H V I E H I V K    Dog
5   Numer               19  H S S F I E I V N I E H V I E H I V K    Dog
6   Numer               19  H S S F I E I V N I E H V I E H I V K    Dog
7   Numer               19  H S S F I E I V N I E H V I E H I V K    Dog
8   Numer               19  H S S F I E I V N I E H V I E H I V K  Human
9   Numer               19  H S S F I E I V N I E H V I E H I V K  Human
10  Numer               19  H S S F I E I V N I E H V I E H I V K  Human
11  Numer               19  H S S F I E I V N I E H V I E H I V K  Human

But I would like that the number of rows for Human is the same for Dog (because, in other words, I would like the same number of data for each class).
Hoping to have been clear, I thank you in advance.
P.S.: My program is an adapted version of this notebook .

Please correct me if I am wrong, you want to print the lines Human 4 times because it is equal to the number of rows that corresponds to Dog, is that right?

Assuming what I said in the comments is true. Here's my attempt, there may be an easier solution but this should do it:
import pandas as pd
import re

# a function to return the count number of animals in the fasta in a dictionary
def count_animals(file_path, separator, animal_index):
  animal_dictionary={}
  with open("Proof.txt") as fasta_file:
    fasta_contents = fasta_file.readlines()
    for line in fasta_contents:
      if line.startswith(">"):
        animal = line.split(separator)[animal_index]
        if animal not in animal_dictionary:
          animal_dictionary[animal] = 1
        else:
          animal_dictionary[animal] = animal_dictionary[animal] + 1
  return animal_dictionary

def read_fasta(file_path, columns) :
  counter = 0  # counter used to print no more than number of times 'Dog'
  animal_d = count_animals(file_path, '|', 1) # execute function
  max_animal =  max(animal_d.values()) # get max number of animals
  min_animal =  min(animal_d.values()) # get min number of animals

  from Bio.SeqIO.FastaIO import SimpleFastaParser
  with open("Proof.txt") as fasta_file :
    records = [] # create empty list
    for title, sequence in SimpleFastaParser(fasta_file): #SimpleFastaParser Iterate over Fasta records as string tuples. For each record a tuple of two strings is returned, the FASTA title line (without the leading ‘>’ character), and the sequence (with any whitespace removed).
      record = []
      title_splits=re.findall(r"[\w']+", title) # Data cleaning is needed

      record.append(title_splits[0])  #First values are ID (Append adds element to a list)
      record.append(len(sequence)) #Second values are sequences lengths
      sequence = " ".join(sequence) #It converts into one line
      record.append(sequence)#Third values are sequences
      # if 'Human' exists in title and the counter is less than min_animal('Dog'), add the record
      if 'Human' in title and counter < min_animal:
        counter += 1 # add +1 to counter

        #Fourth column will contain the species
        record.append("Human")

        records.append(record)
      elif 'Dog' in title:
        record.append("Dog")

        records.append(record)

    # print(records)
  return pd.DataFrame(records, columns = columns) #We have created a function that returns a dataframe

#Now let's use this function by inserting in the first argument the file name (or file path if your working directory is different from where the fasta file is)
#And in the second one the names of columns
data = read_fasta("Proof.txt", columns=["id","sequence_length", "sequence", "animal"])

print(data)

Should give:
      id  sequence_length                               sequence animal
0  Numer               19  H S S F I E I V N I E H V I E H I V K  Human
1  Numer               19  H S S F I E I V N I E H V I E H I V K  Human
2  Numer               19  H S S F I E I V N I E H V I E H I V K  Human
3  Numer               19  H S S F I E I V N I E H V I E H I V K  Human
4  Numer               19  H S S F I E I V N I E H V I E H I V K    Dog
5  Numer               19  H S S F I E I V N I E H V I E H I V K    Dog
6  Numer               19  H S S F I E I V N I E H V I E H I V K    Dog
7  Numer               19  H S S F I E I V N I E H V I E H I V K    Dog

Thank you very much !

An example of contructing a pandas dataframe for fasta on Stack Exchange is here . This is clean code deploying zip outside Biopython. Constructing a pandas dataframe directly via dictionary.
Are you intending to use split method directly on the sequence, such that each residue is a separate column in a dataframe?
Anyway,
title_splits=re.findall(r"[\w']+", title) # Data cleaning is needed

This is a bad idea because you will remove the Genbank unique identifier (its an integer) and this should not leave the sequence.
You will need to switch out of the Fasta object. Your approacha and @user324810's use of
sequence = " ".join(sequence)

is okayish. However, I think it should be '', or "" (no space character). I spotted manually removing fasta, well the method should do this automatically. I am not familiar with Fastaparser, but AlignIO and SeqIO will both automatically process fasta without the need deal with ">". Biopython will have a length method and whilst you have the Biopython object it is preferable to use Biopython (OOP) methods.
from Bio.SeqIO.FastaIO import SimpleFastaParser

This should be at the top of the code.
Overall, you manually open the data, pass it through the Bioparser then immediately dump the object. I do agree to write a pandas dataframe from a Biopython object these conversion to string, but Biopython will provide a method to print the object (pass it to a string) without all the parsing. Some of the parsing is removing key data (Genbank ID).
Personally I would look at the Stackoverflow example and base your code around this.
You clearly have an understanding of pandas, but I get the impression Biopython is something your not familiar with. My suspicion is that Biopython would infact answer the entirety of your question (whatever it might be) without passing to pandas.
pandas is cool and powerful for non-molecular data, but unless your using machine/deep learning (where a dataframe is necessary input) I think Biopython will provide most of the answers.
Important Terminology
The class is FastaIO and the method is SimpleFastaParser. Humans and dogs are instances of this method and are therefore object(s).
Please don't write 'class' of dogs and humans in a formal report because it is incorrect OOP and will cause alarm if your reviewer is a coder. Your code is functional programming rather than OOP. You could say the human/dog fasta object, or just refer to it as human/dog dataframes and avoid any OOP terminology.

Thank you for the suggestions @M__ but SimpleFastaParser does not remove Genbank ID... or at least the FASTA file that I have is with the header like this: >Numer|Animal|OtherInformations . And once SimpleFastaParser, for each record, has returned a tuple of two strings (and the header that is without ">"), re.findall has splitted the header string according to all characters (or at least according to the "|") that I am interested in, I can use this splitting to identify the name of the species and append it to the empty list...

The only thing is in fact that re.findall splits also according to point as separator and I noticed that GenBank IDs are numbers like NP288364.1 ... So re.findall will split NP288364.1 into NP288364 and 1 . Is it so terrible ? I found that the ".1" indicates the "version" of the sequence. Do you know a way to change re.findall so that I can split without considering the points ?

Yes Biopython is new for me but it does not delete the ID... maybe the fault lies with re.findall, isn't it ?

Yes, you should use different parsing such as re.sub explicitly stating "|" such as using a capture "(.+)|.*" and $1, or \1 in re.sub or using or re.final ([\w|\d], string).

Just to add, usually the way to move a sequence object to a neat, decluttered string in Biopython is not to do this manually as done here:
title_splits=re.findall(r"[\w']+", title)
sequence = " ".join(sequence)

both being from SimpleFastaParser(fastafile)
It should be possible and, if so, better to state,
 str(title)
 str(sequence)

If you can issue a print command and the string appears nice and neat on the screen str() will work.
This becomes important for complex sequence/alignment objects, such as SeqIO and AlignIO. The reason it works is __str__ has been declared or available through OOP inheritance.

Thank you, however in fact I realize only now that sequence = " ".join(sequence) is not necessary since SimpleFastaParser already joins sequences into one line and as regards the header I tried with str as you say but it does not work because it leaves it as it is Numer|Human|Hearth and I do not know how to say "append Number". I could write like  title_str=str(title) record.append(title_str[0:6]) but the length of the id number varies for each header...

Fine good to know. str isn't implemented for this given method (they haven't written one). It is a standard method in most OOP, but not all


text:
How many protein sequences do we usually have for each species?
My question is exactly the one in the title. I add a particular example for more clarity. If I want to study the sequence of protein IkBα, do I find only one sequence or is there a database that contains many sequences from many human host species ?
EDIT: The analysis that I am going to do is to use NLP models to apply it to protein sequences to extrapolate, for instance, protein function knowledge.
Here is a reference of a research paper Learning the protein language: Evolution, structure, and function that considers this application by taking into account the problem of low number of sequences data.
So I think I should have at least 10 protein sequences.
Thank you in advance.

@M__ of course, I am going to read Maximillian's response and to answer !

Just to add to Maximilian Press's answer. I happen to know about human immunology and I think Maximilian Press has a good idea.
Genetic studies and immunological studies per se in cellular immunology  tend to be restricted within humans. There are limited exceptions to this namely chimpanzees and humanised mammals, most notably humanised mice. In cellular immunology these mammals retain immunologically distinct behaviour but are sufficiently similar to be useful models.
The protein in question is a B-cell repressor, thats more complicated ...
However, using trans-species homology of the protein of interest for

chimpanzees
humanised mice
regular mice

.. represents a good diversity of proteins to use what appears to be a deep-learning function prediction if this was a traditional cellular immunological study. Any immunologist looking at the biological rationale would recognise the model systems deployed.
Method: Blast Generating diversity for given protein is easy, except for humans. Just use online Blast and restrict the search to the species of interest and 20000 hits (?I forget the upper limit) is the maximum and select for more diverse genes. This would work for chimps and (possibly) mice. Thus blast a human gene and restrict the output to the species of interest. Human gene diversity - this much bigger than an online blast could be targeted in accordance with human genetic diversity. The cellular immune system for humans does not conform to background diversity in general. Some geographic human populations do have a HLA predominance but there are no absolute association.
The final option is local blasting the human database for that gene, that is doable and objectively used to derive representative diversity - for this purpose is complicated and would advise against it.

thanks @M__ . Even if you define your answer as an "addition" to the first one, it was what I was looking for (a general explanation about the database).

Thanks for your feedback @HelpNeederStudent

I'm going to start by assuming that you don't actually want 10 protein sequences from the same species, even though you write this in your comment on M__'s (deleted) answer. This will not be a very informative sample of variation in a given protein, except for extreme examples like antibodies.
To get a good sample of variation you almost certainly want to look across species. The paper that you link to does this by using the UniRef database, they then cluster proteins at the 90% similarity threshold. That is well beyond the variation that you'd expect within a species. See their section "Training datasets":

We train our masked language models on a large corpus of protein sequences, UniRef90 (Suzek et al., 2007), retrieved in July 2018. This dataset contains 76,215,872 protein sequences filtered to 90% sequence identity.

You could use an approach like that in the paper to get similar-ish clusters of proteins whose variation you can then analyze. You can use this, but I am somewhat suspicious of that, without reading the paper in too much depth, because we know rather a lot about orthology/paralogy and how it affects function in these cases. Especially in these very elaborate not very introspectable models I worry about what they will do with inadequately curated data.
So instead I'd probably want to first identify the orthology groups (gene families) of interest. These are generally manually curated groupings of proteins.
Once you have the orthology groups that you care about, you can use a database such as KEGG or OMA to identify the specific orthologs (genes) that belong to the group. Then you will have your dataset.

Thank you Maximilian for your explanation about the paper. However, if one wants to search for motifs, the dataset should be 10 protein sequences from the same species, right ? Even though, as @M___ has explained, "the cellular immune system for humans does not conform to background diversity in general. Some geographic human populations do have a HLA predominance but there are no absolute association".

@HelpNeederStudent I am not sure why you would need them to be from the same species. You will have to clarify that. The whole idea of a gene family or orthology group is that it's a bunch of curated genes/proteins that have the same essential biological function. The people in your paper certainly aren't concerned with species of origin! You could still theoretically do it within a species, I'm just skeptical that you'll find a level of variation that you can do any work with (my professional though non-expert opinion).

Thank you very much !


text:
How to identify genes from a genome assembly of C. Elegans?
I have two full genome assemblies for C. Elegans samples collected from two different geographical areas that I found on WormBase. These are in fasta format. I want to go gene-by-gene and compare the nucleotide sequences corresponding to the gene between the two samples. My goal with this is to count how many single-base differences there are in each gene. However, I'm not sure how to match the sequences in the assembly to known reference sequences (by the way is there a good resource to download all of these at once?) for C. Elegans genes.
My thought was to run BLAST between the reference sequence and each genome and find the match with the best score, if any. Is this a reasonable approach, or is there some better way to do this?

You don't want BLAST here, or at least not regular BLASTn. A better method would be to use tBLASTn to map the known C. Elegans proteins to the translated genome. However, I would recommend using a slightly more sophisticated approach:

Collect fasta sequences of all known C. Elegans proteins (amino acid sequences, not nucleotide). This should be easy enough to do on WormBase. For example, the protein sequences for C. elegans strain VC2010, assembly PRJEB28388, can be found here: ftp://ftp.wormbase.org/pub/wormbase/releases/current-production-release/species/c_elegans/PRJEB28388/c_elegans.PRJEB28388.WS279.protein.fa.gz

Use a tool that can model introns and exons to map a protein into a genome. I haven't worked in this field in more than 10 years now, so it is very likely that there are other programs around today, but back in the day, I would do this using either exonerate or genewise. With exonerate, the command would be:
exonerate -m protein2genome -n 1 -t assembly.fasta -q allproteins.pep > out

That will give you a list of best-hits for each of your input proteins and should accurately map the protein to the genome. Once you have that, you can extract the matches and start comparing.


text:
strategy for generating a genotype likelihood file from large bam files and many samples
I am currently trying to run angsd v0.940 on a whole genome shotgun sequence dataset (122 samples; avg file size: 13G, min: 4G, max: 34G) but it keeps running out memory no matter how I try to tune the parameters. For instance, I tried to change the number of threads from very high values like 75 threads to 2 or 1 following suggestions here: https://github.com/ANGSD/angsd/issues/243). My reference genome size is 1GB and 58,006 contigs. I'm trying to generate a genotype likelihood dataset but I can't think of any other way. Any thoughts would be appreciated!
My command:
angsd -bam sample_bam_list -GL 1 -doMajorMinor 1 -doMaf 1 -doGlf 2 -minMaf 0.05 -SNP_pval 0.01 -minInd 122 -minMapQ 20 -minQ 20 -nThreads 1 -out cbth

Any suggestions would be greatly appreciated. Maybe there is a way to break up the bam files into smaller sizes. Alternatively, maybe generate genotype liklihood files in small batches then merge? It does run successfully on 30 of the 122 samples.

What machine are you running it on? Do you have access to a AWS or GCP instance?

Alternatively, maybe generate genotype likelihood files in small batches then merge? It does run successfully on 30 of the 122 samples.

Yes, that would be my recommendation if you're having trouble with memory. The standard strategy recommended in GATK is similar to this, creating genotype likelihoods individually, then merging for genotyping:
https://github.com/broadinstitute/gatk-docs/blob/master/gatk3-methods-and-algorithms/Best_Practices_for_Variant_Discovery_in_DNAseq.md


text:
Determining what RNAseq data is filtered on volcano plot
I am using RNA seq data to analyze genes via a volcano plot comparing differential gene expression of bacteria with and without antibiotic in R. After having created my plot, I am unsure why some of my values which read "0" TPM in one condition and a TPM value that's not "0" in the other condition were not determined to be differentially expressed. Some of the genes on my volcano plot have this difference in TPM and show up as significant on my plot while others with this "0" value difference are not considered to be significant according to my plot.
Sample of my data:
(UO1_D712_##### represents the locus number, the top column represents the different replicates where "1" is untreated and "2" is treated with antibiotic, the letters next to the numbers represent a different technical replicate):
structure(list(`1` = c("U01_D712_00001",

"U01_D712_00002", "U01_D712_00003",
"U01_D712_00004", "U01_D712_00005", "U01_D712_00006", "U01_D712_00007",
"U01_D712_00008", "U01_D712_00009", "U01_D712_00010", "U01_D712_00011",
"U01_D712_00012", "U01_D712_00013", "U01_D712_00014", "U01_D712_00015",
"U01_D712_00016", "U01_D712_00017", "U01_D712_00018", "U01_D712_00019",
"U01_D712_00020", "U01_D712_00021", "U01_D712_00022", "U01_D712_00023",
"U01_D712_00024", "U01_D712_00025", "U01_D712_00026", "U01_D712_00027",
"U01_D712_00028", "U01_D712_00029", "U01_D712_00030"), `1a` = c("6.590456502",
"6.741694688", "7.110342585", "6.299482433", "2.77173982", "2.470330508",
"3.827125186", "6.267229842", "5.524708663", "2.657913228", "2.87209859",
"7.479820548", "4.980185572", "4.210955199", "0", "4.492492822",
"3.611091371", "3.714270433", "7.455036914", "7.045203025", "6.061860857",
"2.925268313", "6.544077039", "5.747318013", "9.97083444", "9.22523089",
"4.20205383", "6.097040679", "2.621192351", "0"), `1b` = c("6.544454427",
"6.601489488", "7.134224619", "5.814043553", "3.280958379", "2.649180803",
"3.860542083", "6.256648363", "5.380427766", "2.581027705", "3.016132165",
"7.405329447", "4.701503289", "4.073814818", "0", "4.304196924",
"3.515977329", "3.843535649", "7.342972625", "6.966606769", "6.122878624",
"3.007522306", "6.495797641", "5.965431621", "9.828050269", "9.219563915",
"4.065778989", "6.105331066", "3.061209408", "0"), `1c` = c("6.608196006",
"6.743010138", "7.102600793", "6.146518601", "3.555184202", "2.364971542",
"4.034053983", "6.158523627", "5.656051812", "2.658660735", "3.054717455",
"7.392164473", "4.950953264", "4.277770477", "0", "4.507936666",
"3.794842979", "3.610794578", "7.471646548", "7.104792624", "6.484767016",
"3.071205184", "6.584425715", "6.20466015", "9.986342122", "9.282943758",
"4.179958213", "6.219551653", "2.984738345", "0"), `1d` = c("6.547155382",
"6.558328892", "6.992501615", "6.449558793", "3.059801464", "2.418800257",
"3.96498952", "6.013208538", "5.279919645", "2.893295085", "1.750510471",
"7.408735671", "4.9425624", "3.804549986", "0", "4.174835979",
"3.806006888", "3.570390524", "7.641137006", "6.976672494", "6.363030106",
"3.083061726", "6.300910093", "6.007490342", "9.926316442", "9.09671588",
"4.320556917", "6.153860107", "2.877230446", "0"), `1e` = c("6.626724417",
"6.577345176", "7.156821278", "6.296873411", "3.618089702", "2.394444986",
"4.129376392", "6.011715246", "5.31197869", "2.00754706", "2.695493528",
"7.538910448", "4.606060035", "3.909472643", "0", "4.346616047",
"3.468681284", "3.338231445", "7.559599613", "7.1527452", "6.232923513",
"3.108624209", "6.535435309", "6.12922864", "10.13108497", "9.310331313",
"3.959568571", "6.182335537", "2.902736258", "0"), `1f` = c("6.419219179",
"6.650459302", "7.319125725", "5.570869357", "2.962845933", "2.55903176",
"4.087597573", "5.995610538", "5.386268651", "2.750800859", "2.416572678",
"7.579148955", "3.952633067", "3.615227674", "0", "4.562838935",
"3.76942104", "3.920096905", "7.935320749", "7.501470652", "6.13700099",
"3.123910608", "6.035971952", "5.706235015", "9.254751395", "8.379630979",
"4.51391973", "5.6890651", "2.43285316", "0"), `1g` = c("6.553221221",
"6.633949852", "7.182305386", "5.769365973", "2.721354972", "1.668390466",
"4.148367057", "6.240883382", "5.458877133", "1.733842637", "2.723803355",
"7.522249899", "4.149567197", "3.780763096", "0", "4.496306813",
"3.645643535", "3.851001768", "7.678552875", "7.283411279", "6.591585956",
"2.879345378", "6.389427003", "5.911222165", "9.851084493", "9.084575304",
"4.272587776", "5.974762147", "2.98852705", "0"), `2a` = c("9.769887737",
"4.550226652", "1.869021464", "4.944848987", "7.9678549", "2.682865013",
"8.495575559", "2.234521659", "3.667196316", "9.180445037", "3.210107621",
"7.21523691", "5.714579923", "5.423986751", "9.118981459", "8.635701597",
"2.742889473", "3.712618983", "8.006057144", "4.999541279", "10.54351774",
"5.880978085", "7.145433526", "6.982416661", "9.339651188", "5.360835327",
"4.699680905", "3.423826225", "6.408271885", "5.038170992"),
    `2b` = c("10.26397519", "1.945664005", "2.086763158", "2.800904763",
    "8.583418657", "1.536094563", "9.32057547", "3.10685839",
    "4.224502319", "10.1252842", "1.811175407", "5.714439316",
    "6.039142559", "3.833361174", "9.757360286", "9.565906731",
    "1.523640473", "2.315033488", "6.312524363", "4.889986456",
    "10.23020108", "4.848727685", "7.533256999", "7.138160378",
    "10.30380331", "3.955469283", "3.167940742", "3.599655687",
    "4.828945262", "3.701043054"), `2c` = c("9.478481216", "2.789289131",
    "3.393949527", "3.754810933", "8.710777154", "1.806170784",
    "9.150005253", "2.612275457", "4.961073313", "9.802701699",
    "2.933183115", "6.532384958", "6.919449225", "4.432699799",
    "9.715063475", "9.265691356", "1.412064593", "3.330131873",
    "6.665979896", "5.158526421", "9.417365584", "4.899531204",
    "8.173459354", "7.271400938", "9.813068613", "4.384622077",
    "3.700645365", "4.457874829", "5.649440022", "3.531010379"
    ), `2d` = c("8.199795497", "2.565711524", "2.202287889",
    "3.856354444", "7.380224849", "2.192476466", "8.14446837",
    "1.144258862", "3.31447122", "8.713146629", "2.697890381",
    "6.304428859", "5.745291803", "4.898396114", "9.173362747",
    "8.339933849", "3.159678152", "4.094587234", "7.608649692",
    "5.280206424", "10.34630403", "5.098585806", "7.262400625",
    "6.150190905", "9.316698845", "5.073027993", "4.695003229",
    "2.485847024", "5.545300465", "4.350571411"), `2e` = c("8.74935033",
    "3.739489484", "0.217205045", "4.413657999", "8.745525588",
    "3.657060086", "9.279834921", "2.898951179", "4.282874018",
    "9.610485827", "3.561102455", "7.228334332", "6.388491443",
    "4.389908652", "8.781086564", "9.178866581", "2.374603596",
    "3.961037408", "7.864369809", "3.654728044", "10.15284858",
    "5.894439123", "7.68020282", "7.12243523", "9.998637438",
    "5.092174395", "4.111530392", "3.776835632", "5.624523213",
    "4.095011377"), `2f` = c("7.648310926", "4.345557215", "1.986576876",
    "4.99426288", "7.087937177", "2.810917253", "7.77706637",
    "2.62822773", "3.581811188", "8.470225989", "3.335757437",
    "7.416094847", "5.208841128", "5.536128034", "8.255571138",
    "7.993319997", "1.9209089", "3.573861828", "7.318814519",
    "5.233804806", "11.05855833", "6.247720809", "6.673407583",
    "6.029960625", "8.806867591", "5.459208493", "4.001428729",
    "3.609936979", "5.876522973", "4.652839671"), `2g` = c("7.555235468",
    "3.892899549", "1.726443458", "5.304546796", "7.039588042",
    "3.027235295", "7.703852207", "1.753183519", "2.909288568",
    "8.385169315", "3.902707541", "7.523315081", "4.978364017",
    "5.49103181", "8.096218606", "7.944822989", "2.352609608",
    "4.155433517", "7.227355741", "5.532668321", "11.24946953",
    "6.159185473", "6.443203375", "5.931761874", "8.7421732",
    "5.502000205", "4.652883503", "3.458323017", "6.566487449",
    "4.89790353")), row.names = c(NA, 30L), class = "data.frame")

Design matrix:
structure(list(SampleID = c("1a ", "1b", "1c",

"1d", "1e", "1f",
"1g", "2a", "2b", "2c", "2d", "2e", "2f", "2g"), Group = c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), Replicate = c("a",
"b", "c", "d", "e", "f", "g", "a", "b", "c", "d", "e", "f", "g"
)), row.names = c(NA, 14L), class = "data.frame")

Here is how I wrote my code to sample the RNA seq data (note ~/Documents/VOLCANO/R10LB_0ugvs0.5ug.csv refers to the the sample of data above and ~/Documents/VOLCANO/R10LB_0vs0.5_designmatrix.csv refers to the above design matrix):
R10LB_0vs0.5 <- read.csv2("~/Documents/VOLCANO/R10LB_0ugvs0.5ug.csv", sep=",", check.names = F)
R10LB_0vs0.5 <- janitor::remove_empty(R10LB_0vs0.5, which = "cols")
R10LB_0vs0.5_design <- read.csv2("~/Documents/VOLCANO/R10LB_0vs0.5_designmatrix.csv", sep=",")

rownames(R10LB_0vs0.5) <- R10LB_0vs0.5$`1`
R10LB_0vs0.5 <- R10LB_0vs0.5[,-1]

#colnames(R10LB_0ugvs0.5ug) <- R10LB_0ugvs0.5ug[1,]
#R10LB_0ugvs0.5ug <- R10LB_0ugvs0.5ug[-1,]

R10LB_0vs0.5 <- data.matrix(R10LB_0vs0.5)

R10LB_0vs0.5 <- DGEList(counts = R10LB_0vs0.5, group = R10LB_0vs0.5_design$Group)
str(R10LB_0vs0.5)

Group <- as.vector(as.character(R10LB_0vs0.5_design$Group))
Replicate <- as.vector((R10LB_0vs0.5_design$Replicate))
R10LB_0vs0.5_designmatrix <- model.matrix(~0+Group+Replicate)

keep <- filterByExpr(R10LB_0vs0.5)
R10LB_0vs0.5_filter <- R10LB_0vs0.5[keep, keep.lib.sizes = FALSE]

R10LB_0vs0.5_Disp <- estimateDisp(R10LB_0vs0.5_filter, R10LB_0vs0.5_designmatrix)

CONTRASTS <- makeContrasts( Group1vs2 = Group1-Group2,
                            levels = R10LB_0vs0.5_designmatrix)

R10LB_0vs0.5fit <- glmFit(R10LB_0vs0.5_Disp, contrast = CONTRASTS[,1])
R10LB_0vs0.5lrt <- glmLRT(R10LB_0vs0.5fit)
R10LB_0vs0.5TT <- topTags(R10LB_0vs0.5lrt, n=nrow(R10LB_0vs0.5_Disp))
write.csv(R10LB_0vs0.5TT, file = "R10LB_0vs0.5_comparison.csv")
saveRDS(R10LB_0vs0.5TT, file = "R10LB_0vs0.5_comparison.RDS")

Volcano plot:
ggplot(data = R10LB_0vs0.5_comparison, aes(x = logFC, y = PValue)) + geom_point()

p4 <- ggplot(data = R10LB_0vs0.5_comparison, aes(x = logFC, y = -log10(PValue))) + geom_point() + theme_minimal()

#Add vertical lines for logFC thresholds, and one horizontal line for the p-value threshold
p5 <- p4 + geom_vline(xintercept = c(-0.6, 0.6), col = "red") + geom_hline(yintercept = -log10(0.05), col = "red")

# add a column of NAs
R10LB_0vs0.5_comparison$diffexpressed <- "NO"
# if logFC > 0.6 and PValue < 0.05, set as "UP"
R10LB_0vs0.5_comparison$diffexpressed[R10LB_0vs0.5_comparison$logFC > 0.6 & R10LB_0vs0.5_comparison$PValue < 0.05] <- "UP"
# if logFC < -0.6 and PValue < 0.05, set as "DOWN"
R10LB_0vs0.5_comparison$diffexpressed[R10LB_0vs0.5_comparison$logFC < -0.6 & R10LB_0vs0.5_comparison$PValue < 0.05] <- "DOWN"

p4 <- ggplot(data=R10LB_0vs0.5_comparison, aes(x=logFC, y=-log10(PValue), col=diffexpressed)) + geom_point() + theme_minimal()

p5 <- p4 + geom_vline(xintercept=c(-0.6, 0.6), col="red") + geom_hline(yintercept=-log10(0.05), col="red")

p6 <- p5 + scale_color_manual(values=c("blue", "black", "red"))
p6 <- p5 + scale_colour_manual(values = mycolors)

R10LB_0vs0.5_comparison$delabel <- NA
R10LB_0vs0.5_comparison$delabel[R10LB_0vs0.5_comparison$diffexpressed != "NO"] <- R10LB_0vs0.5_comparison$gene_symbol[R10LB_0vs0.5_comparison$diffexpressed != "NO"]

ggplot(data=R10LB_0vs0.5_comparison, aes(x=logFC, y=-log10(PValue), col=diffexpressed, label=delabel)) + geom_point() + theme_minimal() +geom_text()

library(ggrepel)
ggplot(data=R10LB_0vs0.5_comparison, aes(x=logFC, y=-log10(PValue), col=diffexpressed, label=delabel)) +geom_point() + theme_minimal() +geom_text_repel() + scale_color_manual(values=c("blue", "black", "red")) + geom_vline(xintercept=c(-0.6, 0.6), col="red") +geom_hline(yintercept=-log10(0.05), col="red")

This is the volcano plot I currently have with all of my data:

Can you explain how some of these functions filtered my data (i.e. filterByExpr, estimateDisp, makeContrasts, glmFit, glmLRT)? Why are some of my data points that change from 0 TPM in one condition to some value in another condition showing up on my plot while others are not? Or should I accept that my data analysis is valid or "correct"?

I heard that low expression genes are often filtered out in a first past due to a high chance of false positives, so they could be dropped from filterByExpr. High variance within a condition might undermine the statistical significance of the mean difference across conditions. Is there a way to look at the variance of my  data, and is there a threshold for variance between my replicates that causes them to be filtered out by filterByExpr or other functions?

It is a long time that I don't use R. However, I had a sinilar problem once. Using not cartesian coordinates, ggplot cut several points out of the plot. If you plan to continue using R, search about this problem. Otherwise, if you can accede to it go with MATLAB that has the best tool for volcano plot, by my opinion. It never fails.

edgeR functions do not do any filtering, i.e., they do not remove any genes from the results or from the plots. The only filtering was that which you did explicitly yourself by:
keep <- filterByExpr(R10LB_0vs0.5)
R10LB_0vs0.5_filter <- R10LB_0vs0.5[keep, keep.lib.sizes = FALSE]

For your data, filterByExpr will keep any gene that is expressed at a minimal level in at least 7 samples.
You could relax that requirement, to allow genes to be expressed in only 5 or 6 samples say, by reducing the large.n argument to filterByExpr.

[disclaimer: I understand DESeq2, not EdgeR]
Most likely, low expression values (as in raw counts, prior to any normalisation such as TPM). Expression that has high dispersion across all samples is less likely to be biologically significant.
It's easier to work out what's going on with an MA plot, with average gene expression on the X axis, and log2 fold change on the Y axis.
I'm used to DESeq2 results, so I'd guess the construction would be something like this (assuming that the 'baseMean' field exists in EdgeR results):
R10LB_0vs0.5lrt %>%
    ggplot() +
    aes(y = logFC, col = padj, x=log(baseMean)) +
    geom_point()

It's possible that EdgeR doesn't have anything like baseMean in its default results. As a proxy for this, you could use the count totals from the input matrix (e.g. R10LB_0vs0.5 %>% rowSums() %>% log2()).

Do you mean that the TPM values are too low? Because I believe the ones that have been filtered out have a similar value to other "0" tpm genes that were significant. Also if I use the average of two conditions, will the genes I pull out be less significant because the MA plot will not account for variance between replicates? Thanks for the answer


text:
How low does an E-Value need to be for a gene to code for the same protein as query?
I'm looking at a bacterial genome and trying to find specific PEGs (Protein Encoding Genes) by BLAST searching. I'm taking the genes found with the lowest E-value in order to psiBLAST them to check it codes for the same protein. I'm not exactly sure of how high the E-Value needs to be for me to conclude that the protein is not coded for by the subject gene.
Does anyone have any advice on how to make this process a little more efficient/scientific?
Thank you.
(By the way, I'm using BLAST with a BLOUSUM62 matrix.)

Just for the sake of clarity, could you please give more detail about what exactly your inputs are for BLAST?
For the first BLAST, I understand it as you having a list of CDS sequences that you're using as individual BLAST queries to a specific bacterial genome, keeping the lowest e-val hit for each query CDS, creating a list of potential PEGs.
For the PSI-BLAST, I understand it as you using the potential PEGs from the first BLAST as query to a list of proteins that you want to match with PEGs.
Let me know if I'm wrong!

I assume that you mean "is homologous to" by "encodes the same protein as". Identical ("same") proteins should be fairly trivial to identify using a tool such as tblastx.
You should probably not use E-value for this. E-value is similar to a p-value, in that it represents some aspect of the data (in this case a sequence and a database) compared to some null hypothesis (meaning in this case statistics on the database). Therefore, it depends both on the specifics of the data and the hypothesis, and there is no universal threshold that will work across all sequences and all databases. (See discussion of the meaning of e-value.)
E-value does not in any way represent the content of the sequence or its relationship to a sequence in the database. You could have E-values close to zero or E-values > 10 for a perfect sequence match or a very loose sequence match, depending on the input sequence and the database.
However, BLAST output does include other information that might be helpful. You could (for example) look at the length of the alignment and the number of mismatches/identities in the alignment. This will give you a pretty good guess at how similar the two sequences are; unless they are very divergent you would expect them to align across most of the length with somewhat high protein identity or similarity (>25% probably?). This process will still not work 100% of the time, but I think it is closer to your intent than any given approach involving E-values.


text:
Automatize classical debugging in a bash script? (normal debugging by hand)
And thanks in advance. I am debugging a text file which is the one the two inputs on the program RGmatch https://bitbucket.org/pfurio/rgmatch/src/master/.
The code works. Is that some records give division by zero.
I am doing it by hand.
But is there a way by making a bash script, to see which records show the problem?
I do not want to touch the python code.
python rgmatch.py -g ../../hg38.ncbiRefSeq.gtf -b ../../test.bed -o myassoSort.txt

I have tried with classical debugging, but my actual question is how to automatize that?
Also there are no empty records. The bed file was created based on refseq.txt from the latest human assembly. And also the GTF file. On the refseq.txt I processed to obtain a bed file.
I have checked if this works on a single record, it does, but not in some specific ones. The error I get is:
  File "rgmatch.py", line 1311, in <module>
    main()
  File "rgmatch.py", line 287, in main
    run(gtf, dhs, outputfile)
  File "rgmatch.py", line 1153, in run
    pctg_area      = (float(region_overlap)/intron_length)*100
ZeroDivisionError: float division by zero

The record I found so far is this one.
chr19    14087850    14090751
I did grep on both GTF files and Refseq Files and in both files there is information about the record. Because the program search for areas between region of interest, If I put inside just the record only in GTF file and bed Files no error occurs.
grep "14090751" refseq.bed
692  NM_138352.3 chr19   -   14087850    14090751    14088416    14090420    6   14087850,14088670,14088879,14089067,14089784,14090084,  14088575,14088782,14088976,14089362,14090084,14090751,  0   SAMD1   cmpl    cmpl    0,2,1,0,0,0,
grep "14090751" hg38.ncbiRefSeq.gtf
chr19   ncbiRefSeq.2021-02-10   exon    14090085    14090751    .   -   .   gene_id "SAMD1"; transcript_id "NM_138352.3"; exon_number "1"; exon_id "NM_138352.3.1"; gene_name "SAMD1";
chr19   ncbiRefSeq.2021-02-10   5UTR    14090421    14090751    .   -   .   gene_id "SAMD1"; transcript_id "NM_138352.3"; exon_number "1"; exon_id "NM_138352.3.1"; gene_name "SAMD1";

Please include the things I asked for over at [linux.se]. We need to see the actual error message. Also, I call it "classical debugging" but I don't think anyone else will understand the term :) So, please [edit] your question, add the specific error message and then show us the commands that you want to automate. Have you checked it it works with just a single record? Does it work with just the first line of the bed file?

Also, if you still need help, feel free to ping me in the site chat room where we can directly debug the issue.

If you want to just start incrementally adding lines from the bed file until you hit the one that reproduces the problem, you can do something like this:
for i in $(seq 20 30 $(wc -l < ../../test.bed)); do
    head -n $i ../../test.bed > slice.bed
    if python rgmatch.py -g ../../hg38.ncbiRefSeq.gtf -b slice.bed -o slice.out
    then
        true
    else
        echo "Failed when running with the $i first lines."
        break
    fi
done

Explanation
seq 20 30 $(wc -l < ../../test.bed) will iterate over all numbers between 20 and the number of lines in your bed file with increments of 30. Next, the head -n $i ../../test.bed > slice.bed will take the first $i lines of the bed file. This way, the loop will run the rgmatch command on the first 20 lines of the file, then the first 50, then the first 80 and so on as long as the command runs successfully. When the rgmatch command fails, it will report the value of $i to help you find the line that caused the problem.

Thank you @terdon I will edit the answer and try this.. If not I will appreciate very much your help. Right now I have a work meeting. But after that, I will do everything you said to me. Thank you !

De res, @EdmondGeraudAguilar! And do feel free to ping me in chat if you need more help.

It works ! I mean, I would that the script makes like a try-catch then store the results.. but I will handle with this. Than you very much

@EdmondGeraudAguilar it's very easy to extend to whatever you need. But this is why I asked you to edit and explain what commands you want to automate. If you want to keep the results of each run, just change the rgmatch output file: if python rgmatch.py -g ../../hg38.ncbiRefSeq.gtf -b slice.bed -o slice.$i.out; then ...

For those who may use this function be aware that some start coordinates of transcripts are missing in the GTF files that is the problem


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Note that the old DSL1 is no longer supported.