Alignment Data

Software and Chunk Options

library(knitr)
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.1     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.1     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

knitr::opts_chunk$set(
  echo = TRUE,         #Display code chunks
  eval = FALSE,         #Evaluate code chunks
  warning = FALSE,     #Hide warnings
  message = FALSE,    # Hide messages
  fig.width = 6,       #Set plot width in inches
  fig.height = 4,      #Set plot height in inches
  fig.align = "center" #Align plots to the center
)

RNA-seq

Look at GTF

head /home/shared/8TB_HDD_01/sr320/github/deep-dive/D-Apul/data/Amil/ncbi_dataset/data/GCF_013753865.1/genomic.gtf 

#gtf-version 2.2
#!genome-build Amil_v2.1
#!genome-build-accession NCBI_Assembly:GCF_013753865.1
#!annotation-source NCBI Acropora millepora Annotation Release 101
NC_058066.1 Gnomon  gene    1962    23221   .   -   .   gene_id "LOC114963522"; transcript_id ""; db_xref "GeneID:114963522"; gbkey "Gene"; gene "LOC114963522"; gene_biotype "lncRNA"; 
NC_058066.1 Gnomon  transcript  1962    23221   .   -   .   gene_id "LOC114963522"; transcript_id "XR_003825913.2"; db_xref "GeneID:114963522"; gbkey "ncRNA"; gene "LOC114963522"; model_evidence "Supporting evidence includes similarity to: 100% coverage of the annotated genomic feature by RNAseq alignments, including 2 samples with support for all annotated introns"; product "uncharacterized LOC114963522"; transcript_biotype "lnc_RNA"; 
NC_058066.1 Gnomon  exon    23085   23221   .   -   .   gene_id "LOC114963522"; transcript_id "XR_003825913.2"; db_xref "GeneID:114963522"; gene "LOC114963522"; model_evidence "Supporting evidence includes similarity to: 100% coverage of the annotated genomic feature by RNAseq alignments, including 2 samples with support for all annotated introns"; product "uncharacterized LOC114963522"; transcript_biotype "lnc_RNA"; exon_number "1"; 
NC_058066.1 Gnomon  exon    21001   21093   .   -   .   gene_id "LOC114963522"; transcript_id "XR_003825913.2"; db_xref "GeneID:114963522"; gene "LOC114963522"; model_evidence "Supporting evidence includes similarity to: 100% coverage of the annotated genomic feature by RNAseq alignments, including 2 samples with support for all annotated introns"; product "uncharacterized LOC114963522"; transcript_biotype "lnc_RNA"; exon_number "2"; 
NC_058066.1 Gnomon  exon    18711   18775   .   -   .   gene_id "LOC114963522"; transcript_id "XR_003825913.2"; db_xref "GeneID:114963522"; gene "LOC114963522"; model_evidence "Supporting evidence includes similarity to: 100% coverage of the annotated genomic feature by RNAseq alignments, including 2 samples with support for all annotated introns"; product "uncharacterized LOC114963522"; transcript_biotype "lnc_RNA"; exon_number "3"; 
NC_058066.1 Gnomon  exon    1962    2119    .   -   .   gene_id "LOC114963522"; transcript_id "XR_003825913.2"; db_xref "GeneID:114963522"; gene "LOC114963522"; model_evidence "Supporting evidence includes similarity to: 100% coverage of the annotated genomic feature by RNAseq alignments, including 2 samples with support for all annotated introns"; product "uncharacterized LOC114963522"; transcript_biotype "lnc_RNA"; exon_number "4"; 

Getting exon coordinates

/home/shared/hisat2-2.2.1/hisat2_extract_exons.py \
/home/shared/8TB_HDD_01/sr320/github/deep-dive/D-Apul/data/Amil/ncbi_dataset/data/GCF_013753865.1/genomic.gtf \
> ../output/06-Apulcra-hisat/m_exon.tab

head ../output/06-Apulcra-hisat/m_exon.tab

NC_058066.1 1961    2118    -
NC_058066.1 15360   15663   +
NC_058066.1 18710   18774   -
NC_058066.1 21000   21092   -
NC_058066.1 22158   22442   +
NC_058066.1 23084   23220   -
NC_058066.1 24770   25066   +
NC_058066.1 26652   26912   +
NC_058066.1 27071   27358   +
NC_058066.1 27658   27951   +

/home/shared/hisat2-2.2.1/hisat2_extract_splice_sites.py \
/home/shared/8TB_HDD_01/sr320/github/deep-dive/D-Apul/data/Amil/ncbi_dataset/data/GCF_013753865.1/genomic.gtf \
> ../output/06-Apulcra-hisat/m_splice_sites.tab

head ../output/06-Apulcra-hisat/m_splice_sites.tab

NC_058066.1 2118    18710   -
NC_058066.1 15663   22158   +
NC_058066.1 18774   21000   -
NC_058066.1 21092   23084   -
NC_058066.1 22442   24770   +
NC_058066.1 25066   26652   +
NC_058066.1 26912   27071   +
NC_058066.1 27358   27658   +
NC_058066.1 27951   28247   +
NC_058066.1 28534   29196   +

Indexing Genome with splice information

Generic Code

"${programs_array[hisat2_build]}" \
"${genome_fasta}" \
"${genome_index_name}" \
--exon "${exons}" \
--ss "${splice_sites}" \
-p "${threads}" \
2> hisat2-build_stats.txt

head /home/shared/8TB_HDD_01/sr320/github/deep-dive/D-Apul/data/Amil/ncbi_dataset/data/GCF_013753865.1/GCF_013753865.1_Amil_v2.1_genomic.fna

>NC_058066.1 Acropora millepora isolate JS-1 chromosome 1, Amil_v2.1, whole genome shotgun sequence
CCATTTGACTGACTGAGCTGTAATTCGGGCAGTTAGCCCGCTGACGTCAACGTTAGAAAATGCAGTGATAACGCACGCAG
AGTCGTAGACATTAAGAGGGAGGTAAAACTGAAGACGCTGGTACTCTAAAGCAAGCAGCAGCTACCAAACAGGAAACATA
AATACAGTGCCGTTTTCCAGCTTTTATCTCTTTCATTTCCTCATTACTTCTAACAGTAATCGTAACCTTATTTCTCAGTG
AGAGTAAACAGGATTTTAAACTTTTATCGTTATCTTAGGGTTGGGACCTTACTCCATGTCAGACTCAGTTAACTGCCTGT
GAGTGTTTGCGCGCGCTACAGCTACTGGACGAGTTAATTTTGCATACTCCATTTTGCTGcgtgctttttttttatgtttt
atttttatagtAATGCAATGAATCATTGGTCTGTTGCAAGATCATGAATTAACAGTCAAAAcaataacatgtttttcatt
CAGTCAGGCGTTTCGAACGACCActctaccaactgagccaaTTCATTAGAAAAAATGTATGAAcacctaaaccctaaccc
ttaccCTAACACCATTAACTACCAGACTGCTATCAAATTGAATGATATATAGCCCTTACGGAACTACGCCTGAAGCTACC
TTTGAATTTCACAATTGCCGCTCCAACGTTTTATTACTTAGCATTCAATCTCATTCAGGCTGCAGCTTTAGTCCTGCTTC

/home/shared/hisat2-2.2.1/hisat2-build \
/home/shared/8TB_HDD_01/sr320/github/deep-dive/D-Apul/data/Amil/ncbi_dataset/data/GCF_013753865.1/GCF_013753865.1_Amil_v2.1_genomic.fna \
../output/06-Apulcra-hisat/GCF_013753865.1_Amil_v2.1 \
--exon ../output/06-Apulcra-hisat/m_exon.tab \
--ss ../output/06-Apulcra-hisat/m_splice_sites.tab \
-p 40 \
/home/shared/8TB_HDD_01/sr320/github/deep-dive/D-Apul/data/Amil/ncbi_dataset/data/GCF_013753865.1/genomic.gtf \
2> ../output/06-Apulcra-hisat/hisat2-build_stats.txt

head -20 ../output/06-Apulcra-hisat/hisat2-build_stats.txt

Settings:
  Output files: "../output/06-Apulcra-hisat/GCF_013753865.1_Amil_v2.1.*.ht2"
  Line rate: 7 (line is 128 bytes)
  Lines per side: 1 (side is 128 bytes)
  Offset rate: 4 (one in 16)
  FTable chars: 10
  Strings: unpacked
  Local offset rate: 3 (one in 8)
  Local fTable chars: 6
  Local sequence length: 57344
  Local sequence overlap between two consecutive indexes: 1024
  Endianness: little
  Actual local endianness: little
  Sanity checking: disabled
  Assertions: disabled
  Random seed: 0
  Sizeofs: void*:8, int:4, long:8, size_t:8
Input files DNA, FASTA:
  /home/shared/8TB_HDD_01/sr320/github/deep-dive/D-Apul/data/Amil/ncbi_dataset/data/GCF_013753865.1/GCF_013753865.1_Amil_v2.1_genomic.fna
Reading reference sizes

Alignment

Generic code

Hisat2 alignments "${programs_array[hisat2]}" \
-x "${genome_index_name}"-1 "${fastq_list_R1}" \
-2 "${fastq_list_R2}"-S "${sample_name}".sam \
2> "${sample_name}"-hisat2_stats.txt

/home/shared/hisat2-2.2.1/hisat2 \
-x ../output/06-Apulcra-hisat/GCF_013753865.1_Amil_v2.1 \
-p 48 \
-1 /home/shared/8TB_HDD_01/sr320/github/deep-dive/D-Apul/data/SRR8601366_1.fastq \
-2 /home/shared/8TB_HDD_01/sr320/github/deep-dive/D-Apul/data/SRR8601366_2.fastq \
-S ../output/06-Apulcra-hisat/SRR8601366_mil.sam \
2>&1 | tee ../output/06-Apulcra-hisat/hisat2_stats.txt

head ../output/06-Apulcra-hisat/SRR8601366_mil.sam

@HD VN:1.0  SO:unsorted
@SQ SN:NC_058066.1  LN:39361238
@SQ SN:NC_058067.1  LN:37394699
@SQ SN:NC_058068.1  LN:19500522
@SQ SN:NC_058069.1  LN:29333776
@SQ SN:NC_058070.1  LN:31085176
@SQ SN:NC_058071.1  LN:19840543
@SQ SN:NC_058072.1  LN:24214911
@SQ SN:NC_058073.1  LN:17823322
@SQ SN:NC_058074.1  LN:19706214

tail ../output/06-Apulcra-hisat/SRR8601366_mil.sam

SRR8601366.9176811  409 NW_025322628.1  3460    0   3S64M9S =   3460    0   GAACTGACCGTTCTGGGTTGCTTTTCCAACCTGGCCACGGATAGAGACTTCAGCATCGCAGATAGGGAGATCCGAC    EEEA/EEAEEA/EEEEEEEE/EEEEEAEE/EEEEEEEEEEEEAEAEEEEEAEEEEEEEAEEEEEEEEEEEAAAAA6    AS:i:-12    ZS:i:-12    XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:64 YT:Z:UP NH:i:5
SRR8601366.9176811  409 NW_025322628.1  7913    0   3S64M9S =   7913    0   GAACTGACCGTTCTGGGTTGCTTTTCCAACCTGGCCACGGATAGAGACTTCAGCATCGCAGATAGGGAGATCCGAC    EEEA/EEAEEA/EEEEEEEE/EEEEEAEE/EEEEEEEEEEEEAEAEEEEEAEEEEEEEAEEEEEEEEEEEAAAAA6    AS:i:-12    ZS:i:-12    XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:64 YT:Z:UP NH:i:5
SRR8601366.9176811  69  NC_058072.1 4564862 0   *   =   4564862 0   TTCTTCCTGAGCGGCGCCGGGTTCATTTCTTCCGTAACCCGCAATGGTGTTTTTACAACCGTAGTTGTAGTACAT AAAAAEEEE/AEEE/EEEE/EEAEEEEEAEEEEEEEEEEEEEEEEEEEAEEEEEEEEA<EEEEEEEEEEEEEEEA YT:Z:UP
SRR8601366.9176830  77  *   0   0   *   *   0   0   NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN ################################### YT:Z:UP YF:Z:NS
SRR8601366.9176830  141 *   0   0   *   *   0   0   GGGGGGGGGGGGGGGGGGGGGGGGGAAAAAGGGGGTGGGGGGGGGGGGGGGGGGAGGTAAAAAAAAAAAGGGGGGG    66A//A////////A//A///////////////////////6EE//666/6////////E/<//EE//A///////    YT:Z:UP
SRR8601366.9176815  83  NW_025322763.1  182426  60  76M =   182365  -137    ATAAAGGAAACATGTCTCACTATGACACAGACGAACCCTATGGTGAAGAAATAGTGGACATAGGAACGACTTTTCA    EEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAAAAA    AS:i:-10    ZS:i:-11    XN:i:0  XM:i:2  XO:i:0  XG:i:0  NM:i:2  MD:Z:9C0A65 YS:i:-11    YT:Z:CP NH:i:1
SRR8601366.9176815  163 NW_025322763.1  182365  60  70M6S   =   182426  137 GTTTGCTTAGCTGAGCACATTTCCGCTATTGAGGACACTAAAAGACCTAACAAACTGTACCATAAAGGAAACATGT    AAAAAEEEEEEEEEEEEEAEEEEEEEEE/EEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEE    AS:i:-11    XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:24A45  YS:i:-10    YT:Z:CP NH:i:1
SRR8601366.9176823  137 NC_058078.1 9201603 0   70M6S   =   9201603 0   ATAAGCATTGTTTCCAGTTTCTCTAGGGACATACAATGGTCCCAAGAGAAAACAAAAACAATCATTATGCGGGGTT    AAAAAEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEAEEE    AS:i:-6 ZS:i:-6 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:70 YT:Z:UP NH:i:2
SRR8601366.9176823  393 NC_058069.1 28068794    0   70M6S   =   28068794    0   ATAAGCATTGTTTCCAGTTTCTCTAGGGACATACAATGGTCCCAAGAGAAAACAAAAACAATCATTATGCGGGGTT    AAAAAEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEAEEE    AS:i:-6 ZS:i:-6 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:70 YT:Z:UP NH:i:2
SRR8601366.9176823  69  NC_058078.1 9201603 0   *   =   9201603 0   TTTTAATAAAACTTGTTATCCTCCCTAATAGGCCTGTTCGGAAAATACCACAATACTCCCTGTCCGCGCGGCCAGC    AAAAAEEAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE6EEEEEEEEEE<//E//////////    YT:Z:UP

This is a single read alignment record in the SAM file format. Each field in the record represents a specific piece of information about the aligned read. Here’s a breakdown of each field in the given SAM record:

1. SRR8601366.9176811: Read name or identifier, a unique string assigned to the read.

2. 409: SAM bitwise FLAG, an integer value that encodes information about the read’s mapping status, strand, and paired-end information. In this case, the value 409 means the read is part of a paired-end sequence and is aligned in the reverse direction.

3. NW_025322628.1: Reference sequence name or chromosome on which the read is aligned.

4. 3460: 1-based leftmost mapping position of the read on the reference sequence.

5. 0: Mapping quality, an integer value that represents the confidence in the alignment. A value of 0 indicates that the mapping quality is unavailable or not provided.

6. 3S64M9S: CIGAR string, which represents the alignment of the read to the reference. In this case, it means 3 soft-clipped bases (3S), 64 matched bases (64M), and 9 soft-clipped bases (9S) at the end.

7. =: Reference name of the mate/next read. The equal sign (=) indicates that the mate read is aligned to the same reference sequence.

8. 3460: 1-based leftmost mapping position of the mate/next read.

9. 0: Observed template length, the inferred size of the whole template or insert. In this case, the value is 0, indicating the information is not available or not meaningful.

10. GAACTGACCGTTCTGGGTTGCTTTTCCAACCTGGCCACGGATAGAGACTTCAGCATCGCAGATAGGGAGATCCGAC: Read sequence, the actual bases of the read.

11. EEEA/EEAEEA/EEEEEEEE/EEEEEAEE/EEEEEEEEEEEEAEAEEEEEAEEEEEEEAEEEEEEEEEEEAAAAA6: Base quality string, Phred-scaled quality scores for each base in the read sequence, encoded as ASCII characters. 12-20: Optional fields, providing additional information about the alignment:

• AS:i:-12: Alignment score.

• ZS:i:-12: Secondary alignment score.

• XN:i:0: Number of ambiguous bases in the reference.

• XM:i:0: Number of mismatches in the alignment.

• XO:i:0: Number of gap openings.

• XG:i:0: Number of gap extensions.

• NM:i:0: Edit distance to the reference.

• MD:Z:64: String for mismatching positions.

• YT:Z:UP: Platform unit.

• NH:i:5: Number of reported alignments that contain the query in the current record.

This SAM record provides a comprehensive representation of the read’s alignment, base qualities, and various other properties that can be used for downstream analyses.

RNA-seq II

cd ../data
curl -O https://gannet.fish.washington.edu/seashell/snaps/HBXI01.fasta

/home/shared/hisat2-2.2.1/hisat2-build \
-f ../data/HBXI01.fasta \
../output/HBXI01.fasta.index

/home/shared/hisat2-2.2.1/hisat2 \
-x ../output/HBXI01.fasta.index \
-p 4 \
-1 /home/shared/8TB_HDD_01/snow_crab/5010/exp01/5010_4_S4_L003_R1_001.fastq.gz \
-2 /home/shared/8TB_HDD_01/snow_crab/5010/exp01/5010_4_S4_L003_R2_001.fastq.gz \
-S ../output/4_HBXI01.sam

head ../output/4_HBXI01.sam
tail ../output/4_HBXI01.sam

@HD VN:1.0  SO:unsorted
@SQ SN:ENA|HBXI01000001|HBXI01000001.1  LN:452
@SQ SN:ENA|HBXI01000002|HBXI01000002.1  LN:311
@SQ SN:ENA|HBXI01000003|HBXI01000003.1  LN:557
@SQ SN:ENA|HBXI01000004|HBXI01000004.1  LN:1101
@SQ SN:ENA|HBXI01000005|HBXI01000005.1  LN:930
@SQ SN:ENA|HBXI01000006|HBXI01000006.1  LN:373
@SQ SN:ENA|HBXI01000007|HBXI01000007.1  LN:351
@SQ SN:ENA|HBXI01000008|HBXI01000008.1  LN:564
@SQ SN:ENA|HBXI01000009|HBXI01000009.1  LN:1289
A01335:111:H2HTFDSX3:3:1335:27597:34538 133 ENA|HBXI01000955|HBXI01000955.1 3529    0   *   =   3529    0   CTCTCGTGTCCTGCCACCACAACACAAGGACACAGTGCTGGGTGGAATTATGTCATGTAAAACCTCACCATTTTTTTTCTCTTTCAACTGATGGTCTTGTAATCAGGAGCAAAAAACTCTTAAGTCTTGTTGGCAATGAATTTTAATATGCAAAAAG   FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:F,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFF:FFFFF   YT:Z:UP
A01335:111:H2HTFDSX3:3:1335:1226:34554  99  ENA|HBXI01008283|HBXI01008283.1 1514    1   157M    =   1612    255 AAAGAGGTCACAGTGTAGAACAGAAACACCCACGGTTCAGGGTTGGAGTCACCGTCAAGATGGCCCTTTCTTTGTCACGCAACACCACCGATGAAGCCAGATAACCAGTACTAAGAAATGTTTAGTCTGTGGGTGGGGAACGAAGCACACTTAACTA   FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFF:F:FFFFFFFFF,FFFFFFFFFFFFF:FF:   AS:i:-5 ZS:i:-5 XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:15C141 YS:i:0  YT:Z:CP NH:i:4
A01335:111:H2HTFDSX3:3:1335:1226:34554  147 ENA|HBXI01008283|HBXI01008283.1 1612    1   157M    =   1514    -255    AGATAACCAGTACTAAGAAATGTTTAGTCTGTGGGTGGGGAACGAAGCACACTTAACTAACTAACATAAAATTCGTGTATGATTACAATGGGTGGTAAAATATACATATAAGTAACTTACAGCATTCATCATCATGAAGTATTCTCTCAATTGTTCC   FFF,FFFFFFFFF:F:FF:FFFFF,:F,FFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFF:FFF:FF:FFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFF:FFF   AS:i:0  ZS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:157    YS:i:-5 YT:Z:CP NH:i:4
A01335:111:H2HTFDSX3:3:1335:1226:34554  355 ENA|HBXI01008284|HBXI01008284.1 1553    1   157M    =   1651    255 AAAGAGGTCACAGTGTAGAACAGAAACACCCACGGTTCAGGGTTGGAGTCACCGTCAAGATGGCCCTTTCTTTGTCACGCAACACCACCGATGAAGCCAGATAACCAGTACTAAGAAATGTTTAGTCTGTGGGTGGGGAACGAAGCACACTTAACTA   FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFF:F:FFFFFFFFF,FFFFFFFFFFFFF:FF:   AS:i:-5 ZS:i:-5 XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:15C141 YS:i:0  YT:Z:CP NH:i:4
A01335:111:H2HTFDSX3:3:1335:1226:34554  403 ENA|HBXI01008284|HBXI01008284.1 1651    1   157M    =   1553    -255    AGATAACCAGTACTAAGAAATGTTTAGTCTGTGGGTGGGGAACGAAGCACACTTAACTAACTAACATAAAATTCGTGTATGATTACAATGGGTGGTAAAATATACATATAAGTAACTTACAGCATTCATCATCATGAAGTATTCTCTCAATTGTTCC   FFF,FFFFFFFFF:F:FF:FFFFF,:F,FFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFF:FFF:FF:FFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFF:FFF   AS:i:0  ZS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:157    YS:i:-5 YT:Z:CP NH:i:4
A01335:111:H2HTFDSX3:3:1335:1226:34554  355 ENA|HBXI01008285|HBXI01008285.1 1436    1   157M    =   1534    255 AAAGAGGTCACAGTGTAGAACAGAAACACCCACGGTTCAGGGTTGGAGTCACCGTCAAGATGGCCCTTTCTTTGTCACGCAACACCACCGATGAAGCCAGATAACCAGTACTAAGAAATGTTTAGTCTGTGGGTGGGGAACGAAGCACACTTAACTA   FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFF:F:FFFFFFFFF,FFFFFFFFFFFFF:FF:   AS:i:-5 ZS:i:-5 XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:15C141 YS:i:0  YT:Z:CP NH:i:4
A01335:111:H2HTFDSX3:3:1335:1226:34554  403 ENA|HBXI01008285|HBXI01008285.1 1534    1   157M    =   1436    -255    AGATAACCAGTACTAAGAAATGTTTAGTCTGTGGGTGGGGAACGAAGCACACTTAACTAACTAACATAAAATTCGTGTATGATTACAATGGGTGGTAAAATATACATATAAGTAACTTACAGCATTCATCATCATGAAGTATTCTCTCAATTGTTCC   FFF,FFFFFFFFF:F:FF:FFFFF,:F,FFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFF:FFF:FF:FFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFF:FFF   AS:i:0  ZS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:157    YS:i:-5 YT:Z:CP NH:i:4
A01335:111:H2HTFDSX3:3:1335:1226:34554  355 ENA|HBXI01008282|HBXI01008282.1 1511    1   157M    =   1609    255 AAAGAGGTCACAGTGTAGAACAGAAACACCCACGGTTCAGGGTTGGAGTCACCGTCAAGATGGCCCTTTCTTTGTCACGCAACACCACCGATGAAGCCAGATAACCAGTACTAAGAAATGTTTAGTCTGTGGGTGGGGAACGAAGCACACTTAACTA   FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFF:F:FFFFFFFFF,FFFFFFFFFFFFF:FF:   AS:i:-5 ZS:i:-5 XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:15C141 YS:i:0  YT:Z:CP NH:i:4
A01335:111:H2HTFDSX3:3:1335:1226:34554  403 ENA|HBXI01008282|HBXI01008282.1 1609    1   157M    =   1511    -255    AGATAACCAGTACTAAGAAATGTTTAGTCTGTGGGTGGGGAACGAAGCACACTTAACTAACTAACATAAAATTCGTGTATGATTACAATGGGTGGTAAAATATACATATAAGTAACTTACAGCATTCATCATCATGAAGTATTCTCTCAATTGTTCC   FFF,FFFFFFFFF:F:FF:FFFFF,:F,FFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFF:FFF:FF:FFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFF:FFF   AS:i:0  ZS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:157    YS:i:-5 YT:Z:CP NH:i:4
A01335:111:H2HTFDSX3:3:1335:3739:34554  77  *   0   0   *   *   0   0   GTTGGAGCTTATTATTGCAAAAACGACGAAAAAGGTCCTTTGGGAGAAGTTAGATCTACCACAGAAAGAAATTTTTGTGGAAGAATGGAGAAGAAAAAGAAAATAAATAAAATAAATAAAATGAATAAGAAAAGAAAATGGAAAGAAAGAAAATGGA   FFFFFFFFFFFFFFFFFFFFFFF

DNA Methylation Alignment

find ${reads_dir}*_1.fq.gz \
| xargs basename -s _R1_val_1.fq.gz | xargs -I{} ${bismark_dir}/bismark \
--path_to_bowtie ${bowtie2_dir} \
-genome ${genome_folder} \
-p 4 \
-score_min L,0,-0.6 \
--non_directional \
-1 ${reads_dir}{}_R1_val_1.fq.gz \
-2 ${reads_dir}{}_R2_val_2.fq.gz

find *.bam | \
xargs basename -s .bam | \
xargs -I{} ${bismark_dir}/deduplicate_bismark \
--bam \
--paired \
{}.bam

# Sort files for methylkit and IGV

find *deduplicated.bam | \
xargs basename -s .bam | \
xargs -I{} ${samtools} \
sort --threads 28 {}.bam \
-o {}.sorted.bam

# Index sorted files for IGV
# The "-@ 16" below specifies number of CPU threads to use.

find *.sorted.bam | \
xargs basename -s .sorted.bam | \
xargs -I{} ${samtools} \
index -@ 28 {}.sorted.bam

/home/shared/samtools-1.12/samtools tview \
../data/19F_R1_val_1_bismark_bt2_pe.deduplicated.sorted.bam \
../data/GCF_002022765.2_C_virginica-3.0_genomic.fa

tview

WGS

cd ../data
curl -O https://owl.fish.washington.edu/nightingales/C_gigas/F143n08_R2_001.fastq.gz
curl -O https://owl.fish.washington.edu/nightingales/C_gigas/F143n08_R1_001.fastq.gz

cd ../data
curl -O https://gannet.fish.washington.edu/panopea/Cg-roslin/cgigas_uk_roslin_v1_genomic-mito.fa
curl -O https://gannet.fish.washington.edu/panopea/Cg-roslin/cgigas_uk_roslin_v1_genomic-mito.fa.fai
curl -O https://gannet.fish.washington.edu/panopea/Cg-roslin/GCF_902806645.1_cgigas_uk_roslin_v1_genomic-mito.gtf

/home/shared/hisat2-2.2.1/hisat2-build \
-f ../data/cgigas_uk_roslin_v1_genomic-mito.fa \
../output/cgigas_uk_roslin_v1_genomic-mito.index

/home/shared/hisat2-2.2.1/hisat2 \
-x ../output/cgigas_uk_roslin_v1_genomic-mito.index \
-p 4 \
-1 ../data/F143n08_R1_001.fastq.gz \
-2 ../data/F143n08_R2_001.fastq.gz \
-S ../output/F143_cgigas.sam

tail -1 ../output/F143_cgigas.sam

NGSNJ-086:949:GW2303262992nd:2:2678:32597:37059 163 NC_047567.1 32432703    60  150M    =   32432971    394 ATGAGACAATTCAATTAAAAAAAGTATTTTCAAGTACGAATGGTTTCATTTAATGTTTTAAAATGATTGATATAAACCTATGTTATTTATAAATCTACAATTTTGAATATTGCTCATCAACAAATCAACAATACATGGTATGTTTTAATT  :FFFFFFFFFFFFFFF:FF:,FFFFFFFFFFF,FFFFF:FFFFFFFFFFFFF::FFFFFFFF:FFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFF:FFFF:F:FFFFFF:FFFFFF,FFFF:FFFF,FFFFFFFFF:FFFFF:F::,  AS:i:-5 ZS:i:-5 XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:17G132 YS:i:-29    YT:Z:CP NH:i:1

# Convert SAM to BAM, using 40 additional threads
/home/shared/samtools-1.12/samtools view -@ 40 -bS \
../output/F143_cgigas.sam > ../output/F143_cgigas.bam

# Sort the BAM file, using 40 additional threads
/home/shared/samtools-1.12/samtools sort -@ 40 \
../output/F143_cgigas.bam -o ../output/F143_cgigas_sorted.bam

# Index the sorted BAM file (multi-threading is not applicable to this operation)
/home/shared/samtools-1.12/samtools index \
../output/F143_cgigas_sorted.bam

Calling SNPs

gatk

/home/shared/gatk-4.2.5.0/gatk CreateSequenceDictionary \
-R ../data/cgigas_uk_roslin_v1_genomic-mito.fa

/home/shared/samtools-1.12/samtools faidx \
../data/cgigas_uk_roslin_v1_genomic-mito.fa

/home/shared/samtools-1.12/samtools view -H ../output/F143_cgigas_sorted.bam | head

/home/shared/samtools-1.12/samtools addreplacerg \
-r "ID:F143" \
-r "SM:F143" \
-o ../output/F143_cgigas_sorted_rg.bam \
../output/F143_cgigas_sorted.bam

/home/shared/samtools-1.12/samtools index ../output/F143_cgigas_sorted_rg.bam

/home/shared/gatk-4.2.5.0/gatk --java-options "-Xmx40G" HaplotypeCaller \
-R ../data/cgigas_uk_roslin_v1_genomic-mito.fa \
-I ../output/F143_cgigas_sorted_rg.bam \
-O ../output/F143_variants.vcf

tail ../output/F143_variants.vcf

NC_047565.1 33427300    .   T   G   96.64   .   AC=1;AF=0.500;AN=2;BaseQRankSum=-1.769;DP=13;ExcessHet=0.0000;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=7.43;ReadPosRankSum=0.916;SOR=0.223  GT:AD:DP:GQ:PL  0/1:9,4:13:99:104,0,337
NC_047565.1 33427334    .   A   G   417.06  .   AC=2;AF=1.00;AN=2;DP=10;ExcessHet=0.0000;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=34.58;SOR=4.804    GT:AD:DP:GQ:PL  1/1:0,10:10:30:431,30,0
NC_047565.1 33427567    .   C   T   210.96  .   AC=2;AF=1.00;AN=2;DP=5;ExcessHet=0.0000;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=26.49;SOR=3.611 GT:AD:DP:GQ:PL  1/1:0,5:5:15:225,15,0
NC_047565.1 33427572    .   A   G   210.96  .   AC=2;AF=1.00;AN=2;DP=5;ExcessHet=0.0000;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=27.93;SOR=3.611 GT:AD:DP:GQ:PL  1/1:0,5:5:15:225,15,0
NC_047565.1 33427585    .   T   C   210.96  .   AC=2;AF=1.00;AN=2;DP=5;ExcessHet=0.0000;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=29.57;SOR=3.611 GT:AD:DP:GQ:PL  1/1:0,5:5:15:225,15,0
NC_047565.1 33427681    .   A   C   35.48   .   AC=2;AF=1.00;AN=2;DP=1;ExcessHet=0.0000;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=60.00;QD=32.18;SOR=1.609    GT:AD:DP:GQ:PL  1/1:0,1:1:3:45,3,0
NC_047565.1 33427691    .   C   T   58.32   .   AC=2;AF=1.00;AN=2;DP=2;ExcessHet=0.0000;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=60.00;QD=29.16;SOR=0.693    GT:AD:DP:GQ:PL  1/1:0,2:2:6:70,6,0
NC_047565.1 33427824    .   T   C   1243.64 .   AC=1;AF=0.500;AN=2;BaseQRankSum=0.506;DP=33;ExcessHet=0.0000;FS=5.828;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=31.86;ReadPosRankSum=1.166;SOR=0.098  GT:AD:DP:GQ:PL  0/1:3,30:33:36:1251,0,36
NC_047565.1 33427832    .   C   A   1243.64 .   AC=1;AF=0.500;AN=2;BaseQRankSum=-1.128;DP=35;ExcessHet=0.0000;FS=5.828;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=33.08;ReadPosRankSum=1.920;SOR=0.098 GT:AD:DP:GQ:PL  0/1:3,30:33:36:1251,0,36
NC_047565.1 33427865    .   T   C   1029.06 .   AC=2;AF=1.00;AN=2;DP=30;ExcessHet=0.0000;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=30.82;SOR=1.688    GT:AD:DP:GQ:PL  1/1:0,29:29:87:1043,87,0

mpileup (follows text)

/home/shared/samtools-1.12/samtools mpileup --no-BAQ \
--fasta-ref ../data/cgigas_uk_roslin_v1_genomic-mito.fa \
../output/F143_cgigas_sorted.bam

/home/shared/bcftools-1.14/bcftools mpileup --threads 40 --no-BAQ \
--fasta-ref ../data/cgigas_uk_roslin_v1_genomic-mito.fa \
../output/F143_cgigas_sorted.bam \
| /home/shared/bcftools-1.14/bcftools call -mv -Oz \
> ../output/F143_mpile.vcf.gz

/home/shared/bcftools-1.14/bcftools mpileup --threads 40 --no-BAQ \
--fasta-ref ../data/cgigas_uk_roslin_v1_genomic-mito.fa \
../output/F143_cgigas_sorted.bam > ../output/F143_mpileup_output.txt

tail ../output/F143_mpileup_output.txt

NC_001276.1 18215   .   A   <*> 0   .   DP=49;I16=2,45,0,0,1773,68061,0,0,2820,169200,0,0,590,10176,0,0;QS=1,0;MQSBZ=0;FS=0;MQ0F=0  PL  0,141,255
NC_001276.1 18216   .   C   <*> 0   .   DP=49;I16=2,45,0,0,1785,68805,0,0,2820,169200,0,0,551,9435,0,0;QS=1,0;MQSBZ=0;FS=0;MQ0F=0   PL  0,141,255
NC_001276.1 18217   .   C   <*> 0   .   DP=49;I16=2,45,0,0,1773,68061,0,0,2820,169200,0,0,512,8772,0,0;QS=1,0;MQSBZ=0;FS=0;MQ0F=0   PL  0,141,255
NC_001276.1 18218   .   A   <*> 0   .   DP=47;I16=2,43,0,0,1661,63331,0,0,2700,162000,0,0,475,8185,0,0;QS=1,0;MQSBZ=0;FS=0;MQ0F=0   PL  0,135,255
NC_001276.1 18219   .   T   <*> 0   .   DP=43;I16=2,39,0,0,1551,59847,0,0,2460,147600,0,0,442,7668,0,0;QS=1,0;MQSBZ=0;FS=0;MQ0F=0   PL  0,123,255
NC_001276.1 18220   .   G   <*> 0   .   DP=40;I16=2,36,0,0,1440,55740,0,0,2280,136800,0,0,412,7214,0,0;QS=1,0;MQSBZ=0;FS=0;MQ0F=0   PL  0,114,255
NC_001276.1 18221   .   A   <*> 0   .   DP=39;I16=2,35,0,0,1415,55115,0,0,2220,133200,0,0,383,6819,0,0;QS=1,0;MQSBZ=0;FS=0;MQ0F=0   PL  0,111,255
NC_001276.1 18222   .   C   <*> 0   .   DP=39;I16=2,35,0,0,1391,53627,0,0,2220,133200,0,0,354,6482,0,0;QS=1,0;MQSBZ=0;FS=0;MQ0F=0   PL  0,111,255
NC_001276.1 18223   .   C   <*> 0   .   DP=34;I16=2,30,0,0,1204,47022,0,0,1920,115200,0,0,330,6198,0,0;QS=1,0;MQSBZ=0;FS=0;MQ0F=0   PL  0,96,255
NC_001276.1 18224   .   A   <*> 0   .   DP=30;I16=2,26,0,0,1082,42794,0,0,1680,100800,0,0,310,5958,0,0;QS=1,0;MQSBZ=0;FS=0;MQ0F=0   PL  0,84,255

cat ../output/F143_mpileup_output.txt \
| /home/shared/bcftools-1.14/bcftools call -mv -Oz \
> ../output/F143_mpile.vcf.gz

zgrep "^##" -v ../output/F143_mpile.vcf.gz | \
awk 'BEGIN{OFS="\t"} {split($8, a, ";"); print $1,$2,$4,$5,$6,a[1],$9,$10}' | tail

NC_001276.1 12378   G   A   225.417 DP=203  GT:PL   1/1:255,255,0
NC_001276.1 12386   G   T   225.417 DP=217  GT:PL   1/1:255,255,0
NC_001276.1 12399   A   G   225.422 DP=244  GT:PL   1/1:255,255,0
NC_001276.1 12650   T   C   228.416 DP=227  GT:PL   1/1:255,255,0
NC_001276.1 12692   C   T   225.417 DP=208  GT:PL   1/1:255,255,0
NC_001276.1 12969   C   T   225.417 DP=234  GT:PL   1/1:255,255,0
NC_001276.1 13635   G   A   225.417 DP=224  GT:PL   1/1:255,255,0
NC_001276.1 16662   T   G   225.417 DP=233  GT:PL   1/1:255,255,0
NC_001276.1 17992   T   C   225.417 DP=224  GT:PL   1/1:255,255,0
NC_001276.1 18044   C   T   225.417 DP=182  GT:PL   1/1:255,255,0

/home/shared/bcftools-1.14/bcftools call \
-v -c --ploidy 2 ../output/F143_mpile.vcf.gz \
> ../output/F143_mpile_calls.vcf

Copy Number Variation

This code is executing a command-line script using the coverageBed tool from the bedtools2 suite. The bedtools2 suite is a powerful set of tools for performing various genomic analyses, particularly for comparing genomic features in the context of the BED, GFF/GTF, and VCF file formats.

Here’s a breakdown of what this specific command is doing:

1. /home/shared/bedtools2/bin/coverageBed: This is the path to the coverageBed tool, which computes the depth of coverage (number of reads mapped) for each genomic feature in a GFF or BED file.

2. -hist: This flag tells the coverageBed tool to output a histogram of coverage for each feature in the GFF file. The histogram reports the coverage at each depth of coverage (from 0 to the max depth observed).

-hist   
Report a histogram of coverage for each feature in A as well as a summary histogram for _all_ features in A.
Output (tab delimited) after each feature in A:
1) depth
2) # bases at depth
3) size of A
4) % of A at depth

1. -a ../data/cgigas_uk_roslin_v1_gene.gff: This specifies the input GFF file that contains the genomic features of interest (e.g., genes, exons, etc.). In this case, it’s a file named cgigas_uk_roslin_v1_gene.gff located in the ../data/ directory.

2. -b ../output/F143_cgigas_sorted_rg.bam: This specifies the input BAM file that contains the aligned sequencing reads. In this case, it’s a file named F143_cgigas_sorted_rg.bam located in the ../output/ directory.

3. > ../output/F143_cov.txt: This redirects the output of the coverageBed tool to a file named F143_cov.txt located in the ../output/ directory. This file will contain the coverage histogram for each feature in the input GFF file.

In summary, this command is calculating the coverage of sequencing reads (from the BAM file) across genomic features (from the GFF file) and generating a histogram of coverage depths for each feature. The output is saved to a file named F143_cov.txt in the ../output/ directory.

/home/shared/bedtools2/bin/coverageBed \
-hist \
-a ../data/cgigas_uk_roslin_v1_gene.gff \
-b ../output/F143_cgigas_sorted_rg.bam \
> ../output/F143_cov.txt

head ../output/F143_cov.txt

NC_047559.1 Gnomon  gene    9839    11386   .   +   .   ID=gene-LOC117693020;Dbxref=GeneID:117693020;Name=LOC117693020;gbkey=Gene;gene=LOC117693020;gene_biotype=lncRNA 208 1   1548    0.0006460
NC_047559.1 Gnomon  gene    9839    11386   .   +   .   ID=gene-LOC117693020;Dbxref=GeneID:117693020;Name=LOC117693020;gbkey=Gene;gene=LOC117693020;gene_biotype=lncRNA 211 1   1548    0.0006460
NC_047559.1 Gnomon  gene    9839    11386   .   +   .   ID=gene-LOC117693020;Dbxref=GeneID:117693020;Name=LOC117693020;gbkey=Gene;gene=LOC117693020;gene_biotype=lncRNA 213 1   1548    0.0006460
NC_047559.1 Gnomon  gene    9839    11386   .   +   .   ID=gene-LOC117693020;Dbxref=GeneID:117693020;Name=LOC117693020;gbkey=Gene;gene=LOC117693020;gene_biotype=lncRNA 217 1   1548    0.0006460
NC_047559.1 Gnomon  gene    9839    11386   .   +   .   ID=gene-LOC117693020;Dbxref=GeneID:117693020;Name=LOC117693020;gbkey=Gene;gene=LOC117693020;gene_biotype=lncRNA 219 1   1548    0.0006460
NC_047559.1 Gnomon  gene    9839    11386   .   +   .   ID=gene-LOC117693020;Dbxref=GeneID:117693020;Name=LOC117693020;gbkey=Gene;gene=LOC117693020;gene_biotype=lncRNA 221 1   1548    0.0006460
NC_047559.1 Gnomon  gene    9839    11386   .   +   .   ID=gene-LOC117693020;Dbxref=GeneID:117693020;Name=LOC117693020;gbkey=Gene;gene=LOC117693020;gene_biotype=lncRNA 224 1   1548    0.0006460
NC_047559.1 Gnomon  gene    9839    11386   .   +   .   ID=gene-LOC117693020;Dbxref=GeneID:117693020;Name=LOC117693020;gbkey=Gene;gene=LOC117693020;gene_biotype=lncRNA 226 1   1548    0.0006460
NC_047559.1 Gnomon  gene    9839    11386   .   +   .   ID=gene-LOC117693020;Dbxref=GeneID:117693020;Name=LOC117693020;gbkey=Gene;gene=LOC117693020;gene_biotype=lncRNA 227 5   1548    0.0032300
NC_047559.1 Gnomon  gene    9839    11386   .   +   .   ID=gene-LOC117693020;Dbxref=GeneID:117693020;Name=LOC117693020;gbkey=Gene;gene=LOC117693020;gene_biotype=lncRNA 228 1   1548    0.0006460

sort -t $'\t' -k10,10 ../output/F143_cov.txt | tail

NC_047561.1 Gnomon  gene    15572416    15582046    .   -   .   ID=gene-LOC105332243;Dbxref=GeneID:105332243;Name=LOC105332243;gbkey=Gene;gene=LOC105332243;gene_biotype=lncRNA 9999    4   9631    0.0004153
NC_047562.1 Gnomon  gene    26560899    26591530    .   +   .   ID=gene-LOC117690251;Dbxref=GeneID:117690251;Name=LOC117690251;gbkey=Gene;gene=LOC117690251;gene_biotype=lncRNA 9999    2   30632   0.0000653
NC_047564.1 Gnomon  gene    35714434    35730048    .   -   .   ID=gene-LOC105323138;Dbxref=GeneID:105323138;Name=LOC105323138;gbkey=Gene;gene=LOC105323138;gene_biotype=protein_coding 9999    1   15615   0.0000640
NC_047564.1 Gnomon  gene    35720209    35724053    .   +   .   ID=gene-LOC117692367;Dbxref=GeneID:117692367;Name=LOC117692367;gbkey=Gene;gene=LOC117692367;gene_biotype=lncRNA 9999    1   3845    0.0002601
NC_047565.1 Gnomon  gene    61225939    61469303    .   -   .   ID=gene-LOC105330065;Dbxref=GeneID:105330065;Name=LOC105330065;gbkey=Gene;gene=LOC105330065;gene_biotype=protein_coding 9999    1   243365  0.0000041
NC_047565.1 Gnomon  gene    61244095    61397746    .   +   .   ID=gene-LOC117681750;Dbxref=GeneID:117681750;Name=LOC117681750;gbkey=Gene;gene=LOC117681750;gene_biotype=protein_coding 9999    1   153652  0.0000065
NC_047565.1 Gnomon  gene    61316460    61347818    .   -   .   ID=gene-LOC105324053;Dbxref=GeneID:105324053;Name=LOC105324053;gbkey=Gene;gene=LOC105324053;gene_biotype=protein_coding 9999    1   31359   0.0000319
NC_047566.1 Gnomon  gene    24611214    24618430    .   +   .   ID=gene-LOC117682217;Dbxref=GeneID:117682217;Name=LOC117682217;gbkey=Gene;gene=LOC117682217;gene_biotype=lncRNA 9999    1   7217    0.0001386
NC_047567.1 Gnomon  gene    19989691    20052092    .   -   .   ID=gene-LOC105343873;Dbxref=GeneID:105343873;Name=LOC105343873;gbkey=Gene;gene=LOC105343873;gene_biotype=protein_coding 9999    2   62402   0.0000321
NC_047567.1 Gnomon  gene    36938411    36941487    .   +   .   ID=gene-LOC117683948;Dbxref=GeneID:117683948;Name=LOC117683948;gbkey=Gene;gene=LOC117683948;gene_biotype=lncRNA 9999    1   3077    0.0003250