Strand-specific lncRNA sequencing and analysis

RNA extraction, quantification and qualification

Two replicate bee samples were prepared for RNA sequencing. Each sample contained 20 nurse bees of Apis cerana collected from two hives in Miyun, Beijing, in August 2019.

1) A. cerana nurse bees without gut were pooled for RNA sequencing. Total RNA was extracted with TRIzol (Thermo Fisher) according to the manufacturer's instruction.

2) RNA degradation and contamination was monitored on 1% agarose gels.

3) RNA purity was checked by measuring OD₂₆₀/OD₂₈₀ with NanoPhotometer spectrophotometer (IMPLEN).

4) RNA concentration was measured with Qubit RNA Assay Kit in Qubit 2.0 Flurometer (Life Technologies).

5) RNA integrity was assessed with the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies).

Library preparation and strand-specific transcriptome sequencing

1) Ribosomal RNA was removed from total RNA with Globin-Zero Gold rRNA removal kit (Epicenter) according to the manufacture's instruction, and rRNA-free residue was cleaned up by ethanol precipitation. 

2) Sequencing libraries were generated using the rRNA-depleted RNA by the NEBNext Ultra Directional RNA Library Prep Kit (NEB) following the manufacturer's recommendations. In order to select cDNA fragments of preferentially 200~300bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter). 

3) The library preparations were sequenced on an Illumina Hiseq XTen platform and 150bp paired-end reads at Novogene, Beijing. 

4) Paired reads were removed if containing more than 10% Ns or more than 50% bases with low quality (Q<=5) in either read. 

Reads mapping to leucokinin receptor (lkr) gene region with HISAT2

## make hisat index, "Acer_ACSNU-2.0.fa", "Acer_ACSNU-2.0.gtf" were from A. cerana reference genome (ACSNU-2.0, GCF_001442555.1) in NCBI

hisat2_extract_splice_sites.py Acer_ACSNU-2.0.gtf > Acer_ACSNU-2.0.ss

hisat2_extract_exons.py Acer_ACSNU-2.0.gtf > Acer_ACSNU-2.0.exon

mkdir Acer_ACSNU-2.0_hisat2_indexes

hisat2-2.1.0/hisat2-build -p 4 \

    --ss Acer_ACSNU-2.0.ss --exon Acer_ACSNU-2.0.exon \

    Acer_ACSNU-2.0.fa Acer_ACSNU-2.0_hisat2_indexes/Acer_ACSNU-2.0

# mapping to genome, use clean data from one sample as example

hisat2-2.1.0/hisat2 -p 16 --rna-strandness RF \

    -S rep1_hisat_strand.sam \

    -x Acer_ACSNU-2.0_hisat2_indexes/Acer_ACSNU-2.0 \

    -1 FRRL190301455-2a_1.all.clean.fq -2 FRRL190301455-2a_2.all.clean.fq

# covert to .bam

samtools sort -@ 4 -o rep1_hisat_strand.bam rep1_hisat_strand.sam

samtools index rep1_hisat_strand.bam

## merge with lkr gene region

samtools view -b -h rep1_hisat_strand.bam \

    "NW_016019231.1:1287898-1304439" > rep1_hisat_strand_lkr.bam

## first read in pair maps to reverse strand, read is mapped in proper pair, read is paired

samtools view -f 83 \

    -b rep1_hisat_strand_lkr.bam > rep1_hisat_strand_lkr.sorted.83.bam

samtools view -f 163 \

    -b rep1_hisat_strand_lkr.bam > rep1_hisat_strand_lkr.sorted.163.bam

## second read in pair maps to reverse strand, read is paired, read mapped to proper pair

samtools view -f 99 \

    -b rep1_hisat_strand_lkr.bam > rep1_hisat_strand_lkr.sorted.99.bam

samtools view -f 147 \

    -b rep1_hisat_strand_lkr.bam > rep1_hisat_strand_lkr.sorted.147.bam

## merge files

samtools merge rep1_hisat_strand_lkr.sorted.83.163.bam \

    rep1_hisat_strand_lkr.sorted.83.bam rep1_hisat_strand_lkr.sorted.163.bam

samtools merge rep1_hisat_strand_lkr.sorted.99.147.bam \

    rep1_hisat_strand_lkr.sorted.99.bam rep1_hisat_strand_lkr.sorted.147.bam

# Create feature files that only contain the features on the forward or reverse strand only. 

cat Acer_ACSNU-2.0.gff | awk '$7 == "+" { print $0 }' \

    > Acer_ACSNU-2.0_forward.gff

cat Acer_ACSNU-2.0.gff | awk '$7 == "-" { print $0 }' \

    > Acer_ACSNU-2.0_reverse.gff

## separate out reads mapping to positive and negative features from the stranded read files to get sense and antisense transcription

bedtools intersect -b Acer_ACSNU-2.0_forward.gff \

    -abam rep1_hisat_strand_lkr.sorted.83.163.bam \

    > rep1_hisat_strand_lkr.sorted.83.163.sense.bam

bedtools intersect -b Acer_ACSNU-2.0_reverse.gff \

    -abam rep1_hisat_strand_lkr.sorted.83.163.bam \

    > rep1_hisat_strand_lkr.sorted.83.163.antisense.bam

bedtools intersect -b Acer_ACSNU-2.0_forward.gff \

    -abam rep1_hisat_strand_lkr.sorted.99.147.bam \

    > rep1_hisat_strand_lkr.sorted.99.147.antisense.bam

bedtools intersect -b Acer_ACSNU-2.0_reverse.gff \

    -abam rep1_hisat_strand_lkr.sorted.99.147.bam \

    > rep1_hisat_strand_lkr.sorted.99.147.sense.bam

samtools index rep1_hisat_strand_lkr.sorted.83.163.sense.bam

samtools index rep1_hisat_strand_lkr.sorted.99.147.antisense.bam

## view the bam files with the software Integrated Genome Browser (IGB).