gene_x 0 like s 580 view s
Tags: pipeline, RNA-seq
prepare input reads and samplesheet
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_15_min_A/EX_15_min_A_1.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_15_min_A/EX_15_min_A_2.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_15_min_B/EX_15_min_B_1.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_15_min_B/EX_15_min_B_2.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_15_min_C/EX_15_min_C_1.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_15_min_C/EX_15_min_C_2.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_1h_A/EX_1h_A_1.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_1h_A/EX_1h_A_2.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_1h_B/EX_1h_B_1.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_1h_B/EX_1h_B_2.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_1h_C/EX_1h_C_1.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_1h_C/EX_1h_C_2.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_2h_A/EX_2h_A_1.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_2h_A/EX_2h_A_2.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_2h_B/EX_2h_B_1.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_2h_B/EX_2h_B_2.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_2h_C/EX_2h_C_1.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_2h_C/EX_2h_C_2.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_4h_A/EX_4h_A_1.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_4h_A/EX_4h_A_2.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_4h_B/EX_4h_B_1.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_4h_B/EX_4h_B_2.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_4h_C/EX_4h_C_1.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_4h_C/EX_4h_C_2.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_6h_A/EX_6h_A_1.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_6h_A/EX_6h_A_2.fq.gz .
ln -s ./F23A430001245-01_SEPzqcuP/soapnuke/clean/EX_6h_B/EX_6h_B_1.fq.gz .
ln -s ./F23A430001245-01_SEPzqcuP/soapnuke/clean/EX_6h_B/EX_6h_B_2.fq.gz .
ln -s ./F23A430001245-01_SEPzqcuP/soapnuke/clean/EX_6h_C/EX_6h_C_1.fq.gz .
ln -s ./F23A430001245-01_SEPzqcuP/soapnuke/clean/EX_6h_C/EX_6h_C_2.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_Day_4_A/EX_Day_4_A_1.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_Day_4_A/EX_Day_4_A_2.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_Day_4_B/EX_Day_4_B_1.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_Day_4_B/EX_Day_4_B_2.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_Day_4_C/EX_Day_4_C_1.fq.gz .
ln -s ./F23A430001245_SEPiovzP/soapnuke/clean/EX_Day_4_C/EX_Day_4_C_2.fq.gz .
sample,fastq_1,fastq_2,strandedness
EX_15_min_A,EX_15_min_A_1.fq.gz,EX_15_min_A_2.fq.gz,auto
EX_15_min_B,EX_15_min_B_1.fq.gz,EX_15_min_B_2.fq.gz,auto
EX_15_min_C,EX_15_min_C_1.fq.gz,EX_15_min_C_2.fq.gz,auto
Notes and Debugs
#-- NOTE1: STAR cannot work simoutenously twice at the same time! --
#-- NOTE2 --
The GTF file might be corrupted!
Stop at line : NC_019234.1 RefSeq exon 1600 1977 . + 0 ID "cds-HXG45_RS00010"; Parent "gene-HXG45_RS00010"; Note "incomplete%3B partial in the middle of a contig%3B missing N-terminus and C-terminus"; Ontology_term "GO:0015074"; end_range "1977,."; gbkey "CDS"; go_process "DNA integration|0015074||IEA"; inference "COORDINATES: similar to AA sequence:RefSeq:WP_010972061.1"; locus_tag "HXG45_RS00010"; partial "true"; product "DDE-type integrase/transposase/recombinase"; pseudo "true"; start_range ".,1600"; transl_table "11"
Error Message: Cannot find gene_id!
#cut -d$'\t' -f2 RP62A.gtf | sort -u
cmsearch
GeneMarkS-2+
Protein Homology --> ProteinHomology
RefSeq
tRNAscan-SE
#cut -d$'\t' -f3 RP62A.gtf | sort -u
gene
pseudogene
CDS
rRNA
tRNA
direct_repeat
exon
ncRNA
binding_site
region
riboswitch
RNase_P_RNA
sequence_feature
SRP_RNA
tmRNA
#-- BUG1: mamba update picard --
#2.18.27-SNAPSHOT --> bioconda::picard=3.0.0
conda install -c conda-forge mamba
mamba update picard
- ca-certificates 2023.7.22 hbcca054_0 conda-forge
+ ca-certificates 2023.11.17 hbcca054_0 conda-forge/linux-64 Cached
- nextflow 20.04.1 hecda079_1 bioconda
+ nextflow 23.10.0 hdfd78af_0 bioconda/noarch Cached
- openjdk 11.0.1 h516909a_1016 conda-forge
+ openjdk 17.0.3 h58dac75_5 conda-forge/linux-64 Cached
- picard 3.0.0 hdfd78af_0 bioconda
+ picard 3.1.1 hdfd78af_0 bioconda/noarch Cached
mamba update rsem
#-- BUG2: mamba update gffread --
cd /mnt/h1/jhuang/DATA/Data_Marc_RNA-seq_Sepidermidis/work/b4/703a2595f234c7866e49d9a033f943
gffread \
RP62A.gff \
--keep-exon-attrs -F -T \
-o RP62A.gtf
- gffread 0.9.12 0 bioconda
+ gffread 0.12.7 hdcf5f25_3 bioconda/linux-64 Cached
#-- NOTE: If you use nf-core/rnaseq for your analysis please cite: --
#* The pipeline
# https://doi.org/10.5281/zenodo.1400710
#* The nf-core framework
# https://doi.org/10.1038/s41587-020-0439-x
#* Software dependencies
# https://github.com/nf-core/rnaseq/blob/master/CITATIONS.md
#-- BUG3: under rnaseq_2021 open R --
#using R host rather than conda R
mamba remove r-base
#-I"/home/jhuang/miniconda3/envs/rnaseq_2021/lib/R/include"
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("SummarizedExperiment")
BiocManager::install("tximport")
#Error in citation("tximeta") : there is no package called ‘tximeta’
# Execution halted
#sudo apt-get update
#sudo apt-get install libcurl4-openssl-dev
#BiocManager::install(version = "3.16")
BiocManager::install(c("AnnotationDbi", "httr", "BiocFileCache", "biomaRt", "GenomicFeatures", "ensembldb", "curl", "AnnotationHub"))
BiocManager::install(c("biomaRt", "GenomicFeatures", "ensembldb", "AnnotationHub", "tximeta"))
BiocManager::install(c("tximeta"))
#conda update -n base -c defaults conda
#conda update --all
#!/usr/bin/env Rscript
/mnt/h1/jhuang/DATA/Data_Marc_RNA-seq_Sepidermidis/rnaseq/bin/salmon_tximport.r \
NULL \
salmon \
salmon.merged
prepare gtf file
# --gtf_extra_attributes [string] By default, the pipeline uses the `gene_name` field to obtain additional gene identifiers from the input GTF file when
running Salmon. [default: gene_name]
# --gtf_group_features [string] Define the attribute type used to group features in the GTF file when running Salmon. [default: gene_id]
# --featurecounts_group_type [string] The attribute type used to group feature types in the GTF file when generating the biotype plot with
featureCounts. [default: gene_biotype]
# --featurecounts_feature_type [string] By default, the pipeline assigns reads based on the 'exon' attribute within the GTF file. [default: exon]
#--gtf_extra_attributes gene
#--gtf_group_features Parent --featurecounts_group_type gene_biotype --featurecounts_feature_type CDS (outdated since [CHANGE1,2,3])
#under hamm
conda activate rnaseq
ln -s /home/jhuang/Tools/nf-core-rnaseq-3.12.0/ rnaseq
# ------------ gff_to_gtf.py, then modify *.gtf file as follows -------------
#Upload the scripts gff_to_gtf.py and modify_gtf.py.
python3 gff_to_gtf.py RP62A.gff3 RP62A.gtf
#[CHANGE1] ID in gene and Parent in CDS --> gene_id; ID in exon --> transcript_id; Name in gene --> gene_name
#[CHANGE2] "Protein Homology" --> "RefSeq"
#[CHANGE3] CDS --> exon
#[CHANGE4?] --gtf_extra_attributes ||||gene|||| refers to ID=gene-SERP_RS00220;Name=noc;gbkey=Gene;||||gene=noc||||; "gene=" exists in both gene and exon. choose gene as gtf_extra_attributes! (see --gtf_extra_attributes gene in the next command)
[1] python3 gff_to_gtf.py RP62A.gff3 RP62A.gtf
[2] python3 modify_gtf.py RP62A.gtf RP62A_m.gtf
[3] sed -i -e 's/gene_id "rna-/gene_id "gene-/g' RP62A_m.gtf
#NOTE that "ProteinHomology exon" are from "ProteinHomology CDS"
[4] MANUALLY checking why two are missing, because the gene has two CDS/exons.
#(rnaseq) jhuang@hamm:/mnt/h1/jhuang/DATA/Data_Marc_RNA-seq_Sepidermidis/results/genome$ grep ">cds" results/genome/genome.transcripts.fa | sed 's/>//g' - | sort > CDS_2.txt
#2415
#(rnaseq) jhuang@hamm:/mnt/h1/jhuang/DATA/Data_Marc_RNA-seq_Sepidermidis/results/genome$ grep ">rna" genome.transcripts.fa | wc -l
#84
#cut -d$'\t' -f9 CDS_1.txt | cut -d';' -f1 | sed 's/ID=//g' - | sort > CDS_1_.txt
#2417
#2268d2267
#< cds-WP_010959142.1
#2401d2399
#< cds-WP_161384733.1
NC_002976.3 RefSeq gene 423814 424930 . + . gene_id "gene-SERP_RS02220"; gene_name "prfB"; gbkey "Gene"; gene "prfB"; gene_biotype "protein_coding"; locus_tag "SERP_RS02220"; old_locus_tag "SE0421%2CSERP0421"
NC_002976.3 ProteinHomology exon 423814 423885 . + 0 transcript_id "cds-WP_010959142.1"; gene_id "gene-SERP_RS02220"; Dbxref "GenBank:WP_010959142.1"; Name "WP_010959142.1"; Note "programmed frameshift"; Ontology_term "GO:0006415,GO:0003747"; exception "ribosomal slippage"; gbkey "CDS"; gene "prfB"; go_function "translation release factor activity|0003747||IEA"; go_process "translational termination|0006415||IEA"; inference "COORDINATES: similar to AA sequence:RefSeq:WP_010959142.1"; locus_tag "SERP_RS02220"; product "peptide chain release factor 2"; protein_id "WP_010959142.1"; transl_table "11"
NC_002976.3 ProteinHomology exon 423887 424930 . + 0 transcript_id "cds-WP_010959142.1"; gene_id "gene-SERP_RS02220"; Dbxref "GenBank:WP_010959142.1"; Name "WP_010959142.1"; Note "programmed frameshift"; Ontology_term "GO:0006415,GO:0003747"; exception "ribosomal slippage"; gbkey "CDS"; gene "prfB"; go_function "translation release factor activity|0003747||IEA"; go_process "translational termination|0006415||IEA"; inference "COORDINATES: similar to AA sequence:RefSeq:WP_010959142.1"; locus_tag "SERP_RS02220"; product "peptide chain release factor 2"; protein_id "WP_010959142.1"; transl_table "11"
NC_002976.3 RefSeq gene 2256196 2256869 . + . gene_id "gene-SERP_RS10980"; gene_name "SERP_RS10980"; gbkey "Gene"; gene_biotype "protein_coding"; locus_tag "SERP_RS10980"; old_locus_tag "SERP2219"
NC_002976.3 ProteinHomology exon 2256196 2256385 . + 0 transcript_id "cds-WP_161384733.1"; gene_id "gene-SERP_RS10980"; Dbxref "GenBank:WP_161384733.1"; Name "WP_161384733.1"; Note "programmed frameshift"; Ontology_term "GO:0004803"; exception "ribosomal slippage"; gbkey "CDS"; go_function "transposase activity|0004803||IEA"; inference "COORDINATES: similar to AA sequence:RefSeq:WP_064305258.1"; locus_tag "SERP_RS10980"; product "IS6-like element IS257 family transposase"; protein_id "WP_161384733.1"; transl_table "11"
NC_002976.3 ProteinHomology exon 2256385 2256869 . + 2 transcript_id "cds-WP_161384733.1"; gene_id "gene-SERP_RS10980"; Dbxref "GenBank:WP_161384733.1"; Name "WP_161384733.1"; Note "programmed frameshift"; Ontology_term "GO:0004803"; exception "ribosomal slippage"; gbkey "CDS"; go_function "transposase activity|0004803||IEA"; inference "COORDINATES: similar to AA sequence:RefSeq:WP_064305258.1"; locus_tag "SERP_RS10980"; product "IS6-like element IS257 family transposase"; protein_id "WP_161384733.1"; transl_table "11"
--> "cds-WP_010959142.1" and "cds-WP_161384733.1" occur twice!
NC_002976.3 RefSeq gene 423814 424930 . + . gene_id "gene-SERP_RS02220"; gene_name "prfB"; gbkey "Gene"; gene "prfB"; gene_biotype "protein_coding"; locus_tag "SERP_RS02220"; old_locus_tag "SE0421%2CSERP0421"
NC_002976.3 ProteinHomology exon 423814 424930 . + 0 transcript_id "cds-WP_010959142.1"; gene_id "gene-SERP_RS02220"; Dbxref "GenBank:WP_010959142.1"; Name "WP_010959142.1"; Note "programmed frameshift"; Ontology_term "GO:0006415,GO:0003747"; exception "ribosomal slippage"; gbkey "CDS"; gene "prfB"; go_function "translation release factor activity|0003747||IEA"; go_process "translational termination|0006415||IEA"; inference "COORDINATES: similar to AA sequence:RefSeq:WP_010959142.1"; locus_tag "SERP_RS02220"; product "peptide chain release factor 2"; protein_id "WP_010959142.1"; transl_table "11"
NC_002976.3 RefSeq gene 2256196 2256869 . + . gene_id "gene-SERP_RS10980"; gene_name "SERP_RS10980"; gbkey "Gene"; gene_biotype "protein_coding"; locus_tag "SERP_RS10980"; old_locus_tag "SERP2219"
NC_002976.3 ProteinHomology exon 2256196 2256869 . + 0 transcript_id "cds-WP_161384733.1"; gene_id "gene-SERP_RS10980"; Dbxref "GenBank:WP_161384733.1"; Name "WP_161384733.1"; Note "programmed frameshift"; Ontology_term "GO:0004803"; exception "ribosomal slippage"; gbkey "CDS"; go_function "transposase activity|0004803||IEA"; inference "COORDINATES: similar to AA sequence:RefSeq:WP_064305258.1"; locus_tag "SERP_RS10980"; product "IS6-like element IS257 family transposase"; protein_id "WP_161384733.1"; transl_table "11"
--> "cds-WP_010959142.1" and "cds-WP_161384733.1" occur once!
ln -s /home/jhuang/Tools/nf-core-rnaseq-3.12.0/ rnaseq
(rnaseq) nextflow run rnaseq/main.nf --input samplesheet.csv --outdir results --fasta RP62A.fasta --gtf RP62A_m.gtf -profile test_full -resume --max_memory 100.GB --max_time 2400.h --save_reference --aligner star_salmon --skip_rseqc --skip_dupradar --skip_preseq --skip_biotype_qc --skip_deseq2_qc --min_mapped_reads 0 --gtf_extra_attributes gene
#grep ">" genome.transcripts.fa | wc -l --> 2499 (2,529= CDSs (2,445) + RNA(84))
#-- Ben project: human + bacterium + plasmid --
.gff3
python3 gff_to_gtf.py CP009367.gff3 CP009367.gtf
python3 modify_gtf.py CP009367.gtf CP009367_m.gtf
sed -i -e 's/gene_id "rna-/gene_id "gene-/g' CP009367_m.gtf
python3 gff_to_gtf.py NC_019234.gff3 NC_019234.gtf
python3 modify_gtf.py NC_019234.gtf NC_019234_m.gtf
sed -i -e 's/gene_id "rna-/gene_id "gene-/g' NC_019234_m.gtf
python3 genbank_to_gtf.py 1585.gb 1585.gtf
python3 modify_gtf.py 1585.gtf 1585_m.gtf
#python3 correct_gtf.py
python3 transform_gtf.py 1585_m.gtf 1585_m_.gtf
sed -i -e 's/transcript_id \"/transcript_id \"tx-/g' 1585_m_.gtf
sed -i -e 's/exon_id \"/exon_id \"exon-/g' 1585_m_.gtf
#BUG: why 7052 >= expected 7014
#-- MANUALLY: gene positions out of boundary --
#Transcript cds-WP_011100764.1 is out of chromosome NC_019234.1's boundary!
NC_019234.1 RefSeq gene 66549 67571 . + . gene_id "gene-HXG45_RS00005"; gene_name "repA"; gbkey "Gene"; gene "repA"; gene_biotype "protein_coding"; locus_tag "HXG45_RS00005"; old_locus_tag "D743_p54"
NC_019234.1 RefSeq exon 66549 67571 . + 0 transcript_id "cds-WP_011100764.1"; gene_id "gene-HXG45_RS00005"; Dbxref "GenBank:WP_011100764.1"; Name "WP_011100764.1"; gbkey "CDS"; gene "repA"; inference "COORDINATES: similar to AA sequence:RefSeq:WP_011100764.1"; locus_tag "HXG45_RS00005"; product "plasmid replication initiator RepA"; protein_id "WP_011100764.1"; transl_table "11"
#-->
NC_019234.1 RefSeq gene 66549 66845 . + . gene_id "gene-HXG45_RS00005_1"; gene_name "repA"; gbkey "Gene"; gene "repA"; gene_biotype "protein_coding"; locus_tag "HXG45_RS00005"; old_locus_tag "D743_p54"
NC_019234.1 RefSeq exon 66549 66845 . + 0 transcript_id "cds-WP_011100764.1_1"; gene_id "gene-HXG45_RS00005_1"; Dbxref "GenBank:WP_011100764.1"; Name "WP_011100764.1"; gbkey "CDS"; gene "repA"; inference "COORDINATES: similar to AA sequence:RefSeq:WP_011100764.1"; locus_tag "HXG45_RS00005"; product "plasmid replication initiator RepA"; protein_id "WP_011100764.1"; transl_table "11"
NC_019234.1 RefSeq gene 1 726 . + . gene_id "gene-HXG45_RS00005_2"; gene_name "repA"; gbkey "Gene"; gene "repA"; gene_biotype "protein_coding"; locus_tag "HXG45_RS00005"; old_locus_tag "D743_p54"
NC_019234.1 RefSeq exon 1 726 . + 0 transcript_id "cds-WP_011100764.1_2"; gene_id "gene-HXG45_RS00005_2"; Dbxref "GenBank:WP_011100764.1"; Name "WP_011100764.1"; gbkey "CDS"; gene "repA"; inference "COORDINATES: similar to AA sequence:RefSeq:WP_011100764.1"; locus_tag "HXG45_RS00005"; product "plasmid replication initiator RepA"; protein_id "WP_011100764.1"; transl_table "11"
(rnaseq) nextflow run rnaseq/main.nf --input samplesheet.csv --outdir results_NC_019234 --fasta NC_019234.fasta --gtf NC_019234_m.gtf -profile test_full -resume --max_memory 100.GB --max_time 2400.h --save_reference --aligner star_salmon --skip_rseqc --skip_dupradar --skip_preseq --skip_biotype_qc --skip_deseq2_qc --min_mapped_reads 0 --gtf_extra_attributes gene
nextflow run rnaseq/main.nf --input samplesheet.csv --outdir results_CP009367 --fasta CP009367.fasta --gtf CP009367_m.gtf -profile test_full -resume --max_memory 100.GB --max_time 2400.h --save_reference --aligner star_salmon --skip_rseqc --skip_dupradar --skip_preseq --skip_biotype_qc --skip_deseq2_qc --min_mapped_reads 0 --gtf_extra_attributes gene
nextflow run rnaseq/main.nf --input samplesheet.csv --outdir results_GRCh38 --genome GRCh38 -profile test_full -resume --max_memory 100.GB --max_time 2400.h --save_align_intermeds --save_unaligned --aligner star_salmon --pseudo_aligner salmon
nextflow run rnaseq/main.nf --input samplesheet.csv --outdir results_1585 --fasta 1585.fasta --gtf 1585_m_.gtf -profile test_full -resume --max_memory 100.GB --max_time 2400.h --save_reference --aligner star_salmon --skip_rseqc --skip_dupradar --skip_preseq --skip_biotype_qc --skip_deseq2_qc --min_mapped_reads 0
#TODO: look if gene-HXG45_RS00010 exist!
# -- check the intermediate results --
#vim salmon.merged.gene_counts.tsv
gene_id gene_name EX_15_min_A
gene-SERP_RS00195 rpmH 298
gene-SERP_RS00200 rnpA 90
gene-SERP_RS00205 mnmE 1733
#vim salmon.merged.transcript_counts.tsv
tx gene_id EX_15_min_A
cds-WP_000240855.1 gene-SERP_RS00195 298
cds-WP_001832522.1 gene-SERP_RS00200 90
cds-WP_001831768.1 gene-SERP_RS00205 1733
gff_to_gtf.py (code)
import sys
def gff3_to_gtf(gff3_file, gtf_file):
with open(gff3_file, 'r') as gff3, open(gtf_file, 'w') as gtf:
for line in gff3:
if line.startswith('#'):
continue # Skip header lines
parts = line.strip().split('\t')
if len(parts) != 9:
continue # Skip invalid lines
# Extract fields from GFF3
seqname, source, feature, start, end, score, strand, frame, attributes = parts
# Convert attributes to GTF format
attr_pairs = [f.strip() for f in attributes.split(';') if f.strip() != '']
attr_str = '; '.join([f'{key} "{value}"' for part in attr_pairs for key, value in [part.split('=')]])
# Write to GTF file
gtf_line = f'{seqname}\t{source}\t{feature}\t{start}\t{end}\t{score}\t{strand}\t{frame}\t{attr_str};\n'
gtf.write(gtf_line)
if __name__ == "__main__":
if len(sys.argv) != 3:
print("Usage: python script.py <input.gff3> <output.gtf>")
sys.exit(1)
gff3_file = sys.argv[1]
gtf_file = sys.argv[2]
gff3_to_gtf(gff3_file, gtf_file)
genbank_to_gtf.py (code)
import argparse
from Bio import SeqIO
def genbank_to_gtf(genbank_file, gtf_file):
with open(genbank_file, "r") as input_handle, open(gtf_file, "w") as output_handle:
for record in SeqIO.parse(input_handle, "genbank"):
for feature in record.features:
if feature.type == "gene" or feature.type == "CDS":
start = feature.location.start.position + 1
end = feature.location.end.position
strand = '+' if feature.location.strand == 1 else '-'
attributes = []
gene_id = None
if "gene_id" in feature.qualifiers:
gene_id = feature.qualifiers["gene_id"][0]
elif "gene" in feature.qualifiers:
gene_id = feature.qualifiers["gene"][0]
elif "locus_tag" in feature.qualifiers:
gene_id = feature.qualifiers["locus_tag"][0]
if gene_id is not None:
attributes.append(f'gene_id "{gene_id}"')
transcript_id = None
if "locus_tag" in feature.qualifiers:
transcript_id = feature.qualifiers["locus_tag"][0]
attributes.append(f'transcript_id "{transcript_id}"')
attribute_str = "; ".join(attributes) + ";"
gtf_line = f"{record.id}\t.\t{feature.type}\t{start}\t{end}\t.\t{strand}\t.\t{attribute_str}\n"
output_handle.write(gtf_line)
def main():
parser = argparse.ArgumentParser(description="Convert GenBank file to GTF format.")
parser.add_argument("genbank_file", help="Input GenBank file")
parser.add_argument("gtf_file", help="Output GTF file")
args = parser.parse_args()
genbank_to_gtf(args.genbank_file, args.gtf_file)
if __name__ == "__main__":
main()
modify_gtf.py (code)
import sys
def modify_gtf(input_file, output_file):
with open(input_file, 'r') as infile, open(output_file, 'w') as outfile:
for line in infile:
if line.startswith('#'):
outfile.write(line) # Copy comment lines
continue
parts = line.strip().split('\t')
if len(parts) != 9:
continue # Skip invalid lines
record_type = parts[2]
attributes = parts[8]
# Change "Protein Homology" to "RefSeq"
if parts[1] == "Protein Homology":
parts[1] = "ProteinHomology"
new_attributes = []
for attr in attributes.split(';'):
if attr.strip():
key_value = attr.strip().split(' ', 1) # Split only on the first space
if len(key_value) == 2:
key, value = key_value
if (record_type == 'gene' or record_type == 'pseudogene') and key == 'ID':
key = 'gene_id'
elif (record_type == 'gene' or record_type == 'pseudogene') and key == 'Name':
key = 'gene_name'
elif (record_type == 'CDS' or record_type == 'tRNA' or record_type == 'rRNA' or record_type == 'SRP_RNA') and key == 'ID':
key = 'transcript_id'
elif (record_type == 'CDS' or record_type == 'tRNA' or record_type == 'rRNA' or record_type == 'SRP_RNA') and key == 'Parent':
key = 'gene_id'
elif (record_type == 'exon') and key == 'ID':
key = 'exon_id'
elif (record_type == 'exon') and key == 'Parent':
key = 'transcript_id'
new_attributes.append(f'gene_id {value}')
new_attributes.append(f'{key} {value}')
else:
new_attributes.append(attr.strip())
# Change record type from CDS to exon
if record_type == "CDS":
parts[2] = "exon"
parts[8] = '; '.join(new_attributes)
outfile.write('\t'.join(parts) + '\n')
if __name__ == "__main__":
if len(sys.argv) != 3:
print("Usage: python modify_gtf.py <input_file> <output_file>")
sys.exit(1)
input_file = sys.argv[1]
output_file = sys.argv[2]
modify_gtf(input_file, output_file)
correct_gtf.py (code)
input_file = '1585_m.gtf' # Replace with your input file path
output_file = '1585_m_.gtf' # Replace with your desired output file path
with open(input_file, 'r') as infile, open(output_file, 'w') as outfile:
for line in infile:
if line.strip() and not line.startswith('#'): # Skip empty or comment lines
parts = line.strip().split('\t')
if parts[2] == 'gene': # Check if the feature type is 'gene'
attributes = parts[8].split(';')
attributes = [attr.strip() for attr in attributes if 'transcript_id' not in attr] # Remove transcript_id for gene
parts[8] = '; '.join(attributes) + ';' if attributes else '.'
line = '\t'.join(parts) + '\n'
outfile.write(line)
print(f"Processed file saved as {output_file}")
transform_gtf.py (code)
import sys
def transform_data(input_data):
output_data = []
for line in input_data:
fields = line.split()
entry_type = fields[2]
attributes = ' '.join(fields[8:])
attributes_dict = dict(item.split(' ', 1) for item in attributes.split('; '))
gene_id = attributes_dict.get('gene_id').replace('"', '')
transcript_id = attributes_dict.get('transcript_id').replace('"', '')
if not gene_id.startswith("SF028"):
gene_name = f'gene_name "{gene_id}"'
gene_id = f'gene_id "{transcript_id}"'
else:
gene_name = ''
if entry_type == "gene":
gene_line = f"{fields[0]}\t{fields[1]}\tgene\t{fields[3]}\t{fields[4]}\t{fields[5]}\t{fields[6]}\t{fields[7]}\t{gene_id}; {gene_name}".strip()
output_data.append(gene_line)
transcript_line = f"{fields[0]}\t{fields[1]}\ttranscript\t{fields[3]}\t{fields[4]}\t{fields[5]}\t{fields[6]}\t{fields[7]}\t{gene_id}; transcript_id \"{transcript_id}\""
output_data.append(transcript_line)
elif entry_type == "exon":
exon_line = f"{fields[0]}\t{fields[1]}\texon\t{fields[3]}\t{fields[4]}\t{fields[5]}\t{fields[6]}\t{fields[7]}\t{gene_id}; transcript_id \"{transcript_id}\"; exon_id \"{transcript_id}\""
output_data.append(exon_line)
return output_data
def main(input_file, output_file):
with open(input_file, 'r') as file:
input_data = file.readlines()
transformed_data = transform_data(input_data)
with open(output_file, 'w') as file:
for line in transformed_data:
file.write(line + '\n')
if __name__ == "__main__":
if len(sys.argv) != 3:
print("Usage: python script.py <input_file> <output_file>")
sys.exit(1)
input_file = sys.argv[1]
output_file = sys.argv[2]
main(input_file, output_file)
点赞本文的读者
还没有人对此文章表态
没有评论
RNA-seq 2024 Ute from raw counts
Preparing a GTF file from GenBank for bacterial RNA-seq analysis, using the example of WA
© 2023 XGenes.com Impressum