Estimate transposable elements copies and proportions

Dear,

I do have a 2x150 mRNA-seq experiment with 3 drosophila replicates :

file            format  type    num_seqs        sum_len  min_len  avg_len  max_len
A1_R1.fastq.gz  FASTQ   DNA   22,721,029  3,408,154,350      150      150      150
A1_R2.fastq.gz  FASTQ   DNA   22,721,029  3,408,154,350      150      150      150
A2_R1.fastq.gz  FASTQ   DNA   24,093,863  3,614,079,450      150      150      150
A2_R2.fastq.gz  FASTQ   DNA   24,093,863  3,614,079,450      150      150      150
A3_R1.fastq.gz  FASTQ   DNA   24,080,758  3,612,113,700      150      150      150
A3_R2.fastq.gz  FASTQ   DNA   24,080,758  3,612,113,700      150      150      150

One goal of this study is to identify the TE copies of a transposon (called "412") potentially exprimed. And if possible, their proportions. For that, I have a multi-fasta file with the 30 known copies we identified. Here is the multi alignment of these copies :

multialignment

-> they are similar at >98% , about ~7.5kb each and only few snvs/indels could difference them. I know a long read sequencing could fix the problem but we now have only short-reads, and that question came up later in our full study.

I tried some things to response it (with some BLAST outputs but too tricky, star/hisat2 mapping and extract the uniquely mapped reads) , but the most robust seems to be using salmon (or kallisto).

I ran salmon directly on the raw reads with :

for sample in A1 A2 A3
do
  salmon quant \
    -i $HOME/Analyse_DM412/mapping/salmon/DM412_all_isoforms.index \
    -l IU \
    -1 $HOME/Analyse_DM412/fastqReads/${sample}_R1.fastq.gz \
    -2 $HOME/Analyse_DM412/fastqReads/${sample}_R2.fastq.gz \
    -p 8 \
    --validateMappings \
    --seqBias --gcBias \
    -z ${sample}_salmon/${sample}_alignment.sam \
    -o ${sample}_salmon
done

I show the quant results only for the 5 first copies :

==> A1_salmon/quant.sf <==
Name    Length  EffectiveLength TPM NumReads
412-1-chrX-ref|NC_004354.4|:10138270-10145710   7441    7830.687    17337.115815    634.456
412-2-chrX-ref|NC_004354.4|:551477-558981   7505    7872.176    41855.810400    1539.839
412-5-chrX-ref|NC_004354.4|:1033407-1040982 7576    7949.727    0.000000    0.000
412-6-chrX-ref|NC_004354.4|:19136607-19143468   6862    6882.823    6041.067240 194.314
412-4-chrX-ref|NC_004354.4|:5524896-5528701 3806    3632.014    98162.348244    1666.162

==> A2_salmon/quant.sf <==
Name    Length  EffectiveLength TPM NumReads
412-1-chrX-ref|NC_004354.4|:10138270-10145710   7441    7869.034    31721.650122    1214.067
412-2-chrX-ref|NC_004354.4|:551477-558981   7505    7982.440    65734.986659    2552.099
412-5-chrX-ref|NC_004354.4|:1033407-1040982 7576    8055.466    0.000000    0.000
412-6-chrX-ref|NC_004354.4|:19136607-19143468   6862    6997.530    5434.429198 184.954
412-4-chrX-ref|NC_004354.4|:5524896-5528701 3806    3543.578    57323.131025    987.955

==> A3_salmon/quant.sf <==
Name    Length  EffectiveLength TPM NumReads
412-1-chrX-ref|NC_004354.4|:10138270-10145710   7441    7795.757    19247.178710    1156.249
412-2-chrX-ref|NC_004354.4|:551477-558981   7505    7836.288    56558.662685    3415.351
412-5-chrX-ref|NC_004354.4|:1033407-1040982 7576    7917.329    0.000000    0.000
412-6-chrX-ref|NC_004354.4|:19136607-19143468   6862    6866.519    10016.912066    530.025
412-4-chrX-ref|NC_004354.4|:5524896-5528701 3806    3616.071    101893.786777   2839.295

-> It is very between replicates. When I do a PCA based on the coverage with plotPCA, they still go together.

If I do filter the secondary/supplementary + <Q20 alignments, I loose all the signal.. Do you think I could rely on the salmon results?

Thanks for the help, Paillou

star transposons salmon

• 802 views • link

updated 6 months ago by samuel.a.odonnell

• written 6 months ago by pailloufat •

If I do filter the secondary/supplementary + <Q20 alignments, I loose all the signal.. Do you think I could rely on the salmon results?

I can't remember exactly how Salmon quantifies the quality score but I believe it isn't a phred score. Q20 is generally a good filter for the raw reads but I am not sure why you are using this value here. Do you lose the signal if you just remove secondary and supplemental alignments? After that You are already using the --validatemappings parameter so the mapping should pass those thresholds, and a post-mapping quality filter is probably not necessary.

Are the TE copies with large gaps functional? And do you need to track individual copies? otherwise, considering the similarity, can you just generate a single or several consensus sequences to simplify the analysis?

• 0 views • link

written 6 months ago by samuel.a.odonnell

Thanks for the answer,

I do use this value just as a "symbolic" value, since in it is a common value in the analysis pipelines. If I only filter the secondary and supplemental alignments, I will also lose the signal. That's why I don't really understand how the --validatemappings works in Salmon (first time I use it).

TE copies with large gaps are functional yes. The TE in question is 412, known in Drosophila with 2 ORFs. I was also considering one or more consensus, with common SNVs/indels to get some clusters (with a custom script, or if you know a tool to achieve it?): I won't be as precise as with single copy sequences, but it could be a start.

In any way, tt the end, long reads sequencing will be the solution.

• 0 views • link

written 6 months ago by pailloufat •

What exactly is the signal that you lose?

Do you care if different version of 412 are expressed or just that this TE in general? I think long-reads would only really help if you want to isolate expression of different versions of 412.

For the consensus sequence, a classic way to do it (and used often when generating TE libraries for better or worse) is using cd-hit-est e.g. cd-hit-est -i 412.fa -o 412.consensus.fa -d 0 -aS 0.8 -c 0.8 -G 0 -g 1 -b 500 following the 80-80-80 rule. This should cluster most of yours, but it could be adjusted to 90-90-90 if you want to see a few variations. I think using a consensus sequence would be just as precise (of course not exactly but it is a proxy for expression regardless) as you average out coverage across the element and the different elements are apparently very similar.

• 0 views • link

written 6 months ago by samuel.a.odonnell

I mean I lose all my alignments if I do remove secondary and supplemental alignments.

For sure, long-reads is the solution. It is in discussion in our team.

Thanks for the cdhit tip. I digged a bit and I think a kmer approach could be suit :

run a tool like jellyfish or UniqueKMER to generate the kmers on my multi fasta
extract the ones specific to one copy, if they are (I hope)
map the reads to these specific kmers
normalize along the number kmers copy mapping

• 0 views • link

written 6 months ago by pailloufat •

If you are losing all your alignments, I think the filtering is just not working properly. In all cases I know of, each read will be given a primary alignment which should remain after your filtering.

Can you show some stats of your alignments? e.g. samtools stats to look at the number of secondary/supp alignments etc.

Then can you show what command you are using for filtering.

• 0 views • link

written 6 months ago by samuel.a.odonnell

0 answers

No answers yet.