異なる近縁種が交雑することにより誕生した倍数体を異質倍数体という。異質倍数体は、母親となる種と父親となる種の両方のゲノムを同時に持つ。両親のゲノムにはオーソログ関係にある遺伝子が存在する。異種間交雑により両親のゲノムが異質倍数体に取り込まれると、それらの遺伝子は重複して存在するようになる。このような異質倍数体が持つ重複遺伝子をホメオログという。異質倍数体の遺伝子発現量解析において、その遺伝子（ホメオログ）は、母親ゲノムから発現されたのか、それとも父親ゲノムから発現されたのかを、考える必要が出てくる。

異質倍数体の細胞から RNA をサンプリングし、その RNA ライブラリーを RNA-Seq によるシーケンシングしたとする。このサンプルに対して、ホメオログの発現量解析を行う場合、シーケンシングしたリードを、母親ゲノム由来か、父親ゲノム由来かに分類する必要がある。このような分類を行ってくれるソフトウェアとして HomeoRoq (Akama et al., 2014) や EAGLE-RC (Kuo et al., 2020) などがある。HomeoRoq では、異質倍数体のサンプルからシーケンシングされたリードを母親ゲノムと父親ゲノムの両方にマッピングし、両者のマッピング結果（ミスマッチ数）を比較して、そのリードが母親ゲノム由来か父親ゲノム由来かを決定している。以下に、HomeoRoq の使い方の例を示す。

なお、現在、HomeoRoq は公開中止となっているため、これを入手したい場合は、論文に書かれている連絡先にお問い合わせするとよいかもしれない。また、HomeoRoq の開発者らは、ホメオログのリードをより正確に分類できる EAGLE-RC の使用を推奨している。

データの準備

異質倍数体の RNA-Seq データとして、公開されている A. kamchatica のデータを利用する。このデータはペアエンドリードで、3 biological replicate が存在する。FASTQ データは DDBJ からダウンロードできる。

wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/DRA001/DRA001091/DRX012202/DRR013381_1.fastq.bz2
wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/DRA001/DRA001091/DRX012202/DRR013381_2.fastq.bz2

wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/DRA001/DRA001091/DRX012203/DRR013382_1.fastq.bz2
wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/DRA001/DRA001091/DRX012203/DRR013382_2.fastq.bz2

wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/DRA001/DRA001091/DRX012204/DRR013383_1.fastq.bz2
wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/DRA001/DRA001091/DRX012204/DRR013383_2.fastq.bz2

bzip で圧縮された FASTQ ファイルを gzip 形式で圧縮し直す。

bzip2 -cd DRR013381_1.fastq.bz2 | gzip > DRR013381_1.fastq.gz
bzip2 -cd DRR013381_2.fastq.bz2 | gzip > DRR013381_2.fastq.gz

bzip2 -cd DRR013382_1.fastq.bz2 | gzip > DRR013382_1.fastq.gz
bzip2 -cd DRR013382_2.fastq.bz2 | gzip > DRR013382_2.fastq.gz

bzip2 -cd DRR013383_1.fastq.bz2 | gzip > DRR013383_1.fastq.gz
bzip2 -cd DRR013383_2.fastq.bz2 | gzip > DRR013383_2.fastq.gz

FASTQ クオリティチェック

Trimmomatic を利用して、FASTQ ファイル中のクオリティの低い塩基をトリムし、クオリティの低いリードや短いリードを除去する。

java -jar trimmomatic-0.36.jar PE -phred33 \
     DRR013381_1.fastq.bz2 DRR013381_2.fastq.gz \
     DRR013381_1.fastq.qc.gz DRR013381_1.fastq.unpaired.gz \
     DRR013381_2.fastq.qc.gz DRR013381_2.fastq.unpaired.gz \
     ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

java -jar trimmomatic-0.36.jar PE -phred33 \
     DRR013382_1.fastq.bz2 DRR013382_2.fastq.gz \
     DRR013382_1.fastq.qc.gz DRR013382_1.fastq.unpaired.gz \
     DRR013382_2.fastq.qc.gz DRR013382_2.fastq.unpaired.gz \
     ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

java -jar trimmomatic-0.36.jar PE -phred33 \
     DRR013383_1.fastq.bz2 DRR013383_2.fastq.gz \
     DRR013383_1.fastq.qc.gz DRR013383_1.fastq.unpaired.gz \
     DRR013383_2.fastq.qc.gz DRR013383_2.fastq.unpaired.gz \
     ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

マッピング

次に、リードを母親（A. halleri）のゲノムと父親（A. lyrata）のゲノムにマッピングする。両者のゲノムとリファレンスは EBI から入手できる。

# halleri
wget http://ftp.ebi.ac.uk/pub/databases/ena/wgs/public/ba/BASO01.fasta.gz
gunzip BASO01.fasta.gz
wget http://seselab.org/homeoroq/ahalleri_v1_0.gff

# lyrata
wget http://ftp.ebi.ac.uk/pub/databases/ena/wgs/public/ba/BASP01.fasta.gz
gunzip BASP01.fasta.gz
wget http://seselab.org/homeoroq/alyrata_v1_0.gff

マッピングは、より多くのミスマッチ数を許容できる STAR を利用する。まず、ゲノムのインデックスを作成する。

mkdir halleri_index
STAR --runMode genomeGenerate \
     --genomeDir ./halleri_index \
     --genomeFastaFiles BASO01.fasta.gz \
     --runThreadN 8 --limitGenomeGenerateRAM 10000000000

mkdir lyrata_index
STAR --runMode genomeGenerate \
     --genomeDir ./lyrata_index \
     --genomeFastaFiles BASP01.fasta.gz \
     --runThreadN 8 --limitGenomeGenerateRAM 10000000000

次に、STAR を利用してマッピングを行う。RNA-Seq データが 3 セットあるので、ここではシェルコマンドの for ループで同じ処理を回す。

seq_libs=(DRR013381 DRR013382 DRR013383)
genomes=(halleri lyrata)

for seq_lib in ${seq_libs[@]}; do
    fq1=${seq_lib}_1.qc.fq.gz
    fq2=${seq_lib}_2.qc.fq.gz

    for genome in ${genomes[@]}; do
        tmp_dir=./tmp/star/${seq_lib}__on__${genome}
        mkdir -p ${tmp_dir}
        cd ${tmp_dir}

        STAR --runThreadN 8 \
             --genomeDir ${genome}_index \
             --readFilesIn ${fq1} ${fq2} \
             --readFilesCommand zcat \
             --genomeLoad NoSharedMemory

        samtools-0.1.19/samtools view -bS Aligned.out.sam > Aligned.out.bam
        samtools-0.1.19/samtools sort Aligned.out.bam ${seq_lib}__on__${genome}
        samtools-0.1.19/samtools index ${seq_lib}__on__${genome}.bam

        cd -
    done
done

ホメオログリードの分類

リードを A. halleri と A. lyrata のゲノムにマッピングしたのち、HomeoRoq を使用して、各リードが A. halleri 由来か、A. lyrata 由来かを分類する。

for seq_lib in ${RNASEQ_LIBS[@]}; do
    map_on_amara=${seq_lib}__on__camara
    map_on_rivularis=${seq_lib}__on__crivularis

    tmp_dir=${TMP_DIR}/readclassify/${seq_id}
    mkdir -p ${tmp_dir}
    cd ${tmp_dir}

    python ${SCRIPTDIR}/rc_filter.py ${BAM_DIR}/${map_on_amara}.bam \
                                     ${tmp_dir}/${map_on_amara}.filtered.bam \
                                     ${tmp_dir}/${map_on_amara}.error.bam

    python ${SCRIPTDIR}/rc_filter.py ${BAM_DIR}/${map_on_rivularis}.bam \
                                     ${tmp_dir}/${map_on_rivularis}.filtered.bam \
                                     ${tmp_dir}/${map_on_rivularis}.error.bam

    python ${SCRIPTDIR}/read_classify.py ${tmp_dir}/${map_on_amara}.filtered.bam \
                                         ${tmp_dir}/${map_on_rivularis}.filtered.bam \
                                         ${tmp_dir}/${map_on_amara}.filtered.rc \
                                         ${tmp_dir}/${map_on_rivularis}.filtered.rc


    for tag in ${TAGS[@]}; do
        fname=${tmp_dir}/${map_on_amara}.filtered.rc_${tag}
        ${BIN}/samtools-0.1.19/samtools sort ${fname}.bam ${fname}.sorted
        ${BIN}/samtools-0.1.19/samtools index ${fname}.sorted.bam
        mv ${fname}.sorted.bam ${RC_DIR}/
        mv ${fname}.sorted.bam.bai ${RC_DIR}/

        fname=${tmp_dir}/${map_on_rivularis}.filtered.rc_${tag}
        ${BIN}/samtools-0.1.19/samtools sort ${fname}.bam ${fname}.sorted
        ${BIN}/samtools-0.1.19/samtools index ${fname}.sorted.bam
        mv ${fname}.sorted.bam ${RC_DIR}/
        mv ${fname}.sorted.bam.bai ${RC_DIR}/
    done
done

これにより、各 RNA-Seq データセットに対して、_orig、_common、および_other などの BAM ファイルが生成される。DRR013381 を例にすると、以下のようなファイルが生成される。

HomeoRoq の分類が終わると、上のように母親ゲノム由来と父親ゲノム由来のリードがそれぞれ別々の BAM ファイルに保存される。このあと、featureCounts あるいは HTSeq などを用いて、カウントデータを取得することができる。