miRNA alignment with Bowtie2

Question

[Chironex]

Hi, I was trying to align my microRNAs using Bowtie2 but I get this result:

bowtie2 --local -p 8 -q --phred33 -D 20 -R 3 -N 0 -L 8 -i S,1,0.50 -x miRNA_mature_index -U ../trimmed_reads/trimming_bbduk/H1_trim.fastq.gz -S H1_mature.sam
9637480 reads; of these:
  9637480 (100.00%) were unpaired; of these:
    9567222 (99.27%) aligned 0 times
    1182 (0.01%) aligned exactly 1 time
    69076 (0.72%) aligned >1 times
0.73% overall alignment rate

and I used for trimming this command:

bbduk.sh -Xmx2g in="$input_file" out="$output_file" literal=GCGCCGCTGGTGTAGTGGTATCATGCAAGATTCCAACTGTAGGCACCATC ktrim=r k=7 &> "${file_name}.trim.log"

If I set k=15, alignment rate is 0%

Any advices?

[GenoMax]

You want ungapped alignments so it may be better to try bowtie v.1.x as well.

What is the size of the sequences left over after you do that trimming? Are they between 20-26 bp? If that is your miRNA adapter then you should choose only those sequences that have that adapter.

[Chironex]

Hi, the sequence length after trimming is this:

 740239 13
 682748 14
 646079 15
 585788 16
 550914 17
 513787 18
 443601 19
 400977 20
 372279 21
 349621 22
 249551 23
 208853 24
 183582 25
 151978 26
 133445 27
 107530 28
  99556 29
  88625 30
  70126 31
  59516 32
  52467 33
  47302 34
  41292 35
  33493 36
  28873 37
  25342 38
  21428 39
  19265 40
  17061 41
  14237 42
  12688 43
  10725 44
   9016 45
   7999 46
   6859 47
   6257 48
   5043 49
   4735 50
   3896 51
   3364 52
   3101 53
   2729 54
   2085 55
   2223 56
   1678 57
   1740 58
   1527 59
   1317 60
   1089 61
    959 62
    748 63
    674 64
    689 65
    586 66
    502 67
    413 68
    441 69
    418 70
    348 71
    310 72
    311 73
    262 74
    218 75
    229 76
    189 77
    157 78
    139 79
    114 80
    135 81
     93 82
     42 83
     64 84
     68 85
     47 86
     36 87
     48 88
     53 89
     28 90
     21 91
     45 92
     43 93
     41 94
   2381 101

I'm not able to find a proper trimming for data. However, if I run bowtie v.1.3, results are better:

bowtie --wrapper basic-0 -n 0 -l 15 -e 99999 -k 200 --best --strata -x miRNA_mature_index ../trimmed_reads/trimming_bbduk/H1_trim.fastq.gz -S H1_mature_bbduk.sam
# reads processed: 9637480
# reads with at least one alignment: 732100 (7.60%)
# reads that failed to align: 8905380 (92.40%)
Reported 18004343 alignments

However I'm afraid from trimming step. if I use cutadapt, cutadapt -u 4 -m 15 -q 20 -e 0.1 -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -g ^CAGTCG -g ^TGACTC -g ^GCTAGA -g ^ATCGAT -o H1_R1.fastq.gz there is 0% alignment rate!

[GenoMax]

You should figure out the exact miRNA adapter used in this dataset by finding which kit was used. Looks for read that contain that adapter. After trimming you should be left with reads that will be small 20-25bp. Those would be ones to use.

[Chironex]

Dear GenoMax , I am asking to the guy lab the kit used. However, I will show you the fastq first lines.

@VH00643:236:AACFYWKHV:1:1101:46597:1076 1:N:0:CACTGAACCG+TCCTAAGGTT
GGTAAAGGCTCACCAAGGCGAACTGTAGGCACCATCAATTAGCCAATAGATAGATCGGAAGAGCACACGTCTGAACTCCAGTCACCACTGAACCGATCTCG
+
CC;---CCCC;CC;;CCCCCC;CCCCCCCCC;C-C-CCCC-C-CCCC-C-C;C-CCC;CCC;CCC;CCCCCCCCCCC;CCCCCC-C-CCC-;;CCC;-;-C
@VH00643:236:AACFYWKHV:1:1101:48680:1151 1:N:0:CACTGAACCG+TCCTAAGGTT
GTAGGAGTCTTGCGCAACTGTAGGCACCATCAATGCCACATAAAAAAGATCGGAAGAGCACACGTCTGAACTCCAGTCACCACTGAACCGATGTGGGTTGG
+
;--CC-C--CCCCCCCCC;CCCCCCC;-CCC-CCC-C;C-C--;;CCC;CCCCCCC-CCC;-;CC-CCCCCCC--CCC;C---CC;-C-CCC--CCC-C-C
@VH00643:236:AACFYWKHV:1:1101:44192:1208 1:N:0:CACTGAACCG+TCCTAAGGTT
GTGGCTAGAGGAACTGTAGGCACCATCAATTCTTGCCAACCTAGATCGGAAGAGCACACGTCTGAACTCCAGTCACCACTGAACCGATCTGGGTTGGGGTG
+
;-CC-C-C-CC;CCCCC;CC;CCC;C-CCCCCCCCCCCCC-CCC;CCCCCCC;CCCCCCCCCCCCCCC;CCCCCCCC;CCC;CC-CCC---C-CC-C-C-C
@VH00643:236:AACFYWKHV:1:1101:48661:1208 1:N:0:CACTGAACCG+TCCTAAGGTT
TGAAGCCAGAGAACTGTAGGCACCATCAATTTCCAGCGCTACAGATCGGAAGAGCACACGTCTGAACTCCAGTCACCACTGAACCGATCTGGGTTTGCGTG
+
-C--C;;-C-C-C;CCC-CCCCCC-CCCCCCCCCCCCC;CC;-CCCCCCCCC;C;-;CCCCCCC;CCCCC-CCCCCC;CCC;;C-C-C---C--C-C-C;C
@VH00643:236:AACFYWKHV:1:1101:42582:1227 1:N:0:CACTGAACCG+TCCTAAGGTT
TTCAGTGGGGTGGCGAACTGTAGGCACCATCAATCCCGATTCCCGAAGATCGGAAGAGCACACGTCTGAACTCCAGTCACCACTGAACCGATCTCGGTTTG
+
-C;;CCCCCCCCCCCC-CCCCCCCCCCC;CCCCC;CCCCCCCC;CCCCCC;;CCCC-C-;;CCCCC;CCCCCCCCCC;CCC-CCCC;;CCCCCC-C--CC;
@VH00643:236:AACFYWKHV:1:1101:48869:1227 1:N:0:CACTGAACCG+TCCTAAGGTT
GGTGGAGTGATTTGTCTGGTTAATTCCGTTAACGAACGAGACCTTAACCTAACTGTAGGCACCATCAATAACATCGAGGTGTGTTCGGAAGTGCACACGTC
+
;C;CC-CCC-CCCCCCCCCCC;-CCCCCCC-C;C;CCCCC;--CCCC--C;;CCCCCCC-;-C-CC-CC--;;CCC--C-C-C-C-CCC-C-C----;CC-
@VH00643:236:AACFYWKHV:1:1101:56973:1227 1:N:0:CACTGAACCG+TCCTAAGGTT
GCTAAGCAGGGTCGGGCCTAACTGTAGGCACCATCAATGACGCATAACTAAGATCGGAAGAGCACACGTCTGAACTCCAGTCACCACTGAACCGATCTGGG
+
;-C--C--CCCCCCCCCCCCCCCCC-CCCCCCCCCCC-C;CCCCCCCCCCCC-C;;;CCC-CCCCCCCCC;CCCCCCCCCCC;CCCC;CC-C;CCC;C-C-
@VH00643:236:AACFYWKHV:1:1101:45479:1246 1:N:0:CACTGAACCG+TCCTAAGGTT
GTGGGAGGGGCCGTAACTGTAGCACCATCAATAGACATTAAGCAAGATCGGAAGAGCACACGTCTGAACTCCAGTCACCACTGAACGGATGTGGGTTGGGG
+
C;CCCCCCCCCCCC-CCCCCCCCCCCCCCCCC-C-CCC;CCCCCCC-CCCCCCCCCCC-C-CCCCCCCCCCCCCCC;-;-CCCCC--C-C--CC-CC-;;C
@VH00643:236:AACFYWKHV:1:1101:56841:1246 1:N:0:CACTGAACCG+TCCTAAGGTT
GGGGCTTGGAAGAATCAGCGGGGAACTGTAGGCACCATCAATAGCACAGCATGAAGATCGGAAGAGCACACGTCTGAACTCCTGTCACCACTGAACGGATG
+
-CCC-;CCC-CC-CCC-CCCCCCCCCCCCCCCC;CCCCCCC-;CCCCCCC-CCCCC-CCCCCCCCC-CCCCCCCCCCCCCCC-CC--;C-CCC;---C-C;
@VH00643:236:AACFYWKHV:1:1101:59719:1284 1:N:0:CACTGAACCG+TCCTAAGGTT
GTGTCCTCCGGGAAGTCCTGAACTGTAGGCACCATCAATGGCCTCTACTCGAGATCGGAAGAGCACACGTCTGAACTCCAGTCACCACTGAACCGATGTGG
+
C-C-;;CC-CCC;CCCCCCCCCCCCC-CCC;CCCCC-CCCCCCCCCCCC;CCC-C;CCCCCCC;;CC;CCCCCCCCCC-CCCCC;C-CCCCCC-CCC---C
@VH00643:236:AACFYWKHV:1:1101:35690:1303 1:N:0:CACTGAACCG+TCCTAAGGTT
TCCAAAGCTCCTCCGCGCGACGTGTTGCGAACTGTAGGCACCATCAATCCTCCCGGTCCTAGATCGGAAGAGCACTCGTCTGAACTCCAGTCACCACTGAA
+

It seems that some sequences contain CACTGAACCG that is present in the header. They did do demultiplexing. Should I remove also this sequence? Please give me some advice, I'm little confused. With bbduk command, it removes completely the illumina adapters (FASTQC report), but read length are between 10-101 bp.

[GenoMax]

You can try to take a selection of reads and see if you are able to figure out what the adapter could be (programs like fastp may be able to do this).

Example

[Chironex]

thank you so much.

Detecting adapter sequence for read1...
>Illumina TruSeq Adapter Read 1
**AGATCGGAAGAGCACACGTCTGAACTCCAGTCA**

Read1 before filtering:
total reads: 17419022
total bases: 1759321222
Q20 bases: 1621149200(92.1463%)
Q30 bases: 1562425305(88.8084%)

Read1 after filtering:
total reads: 17411386
total bases: 804356780
Q20 bases: 787222714(97.8698%)
Q30 bases: 766504729(95.2941%)

Filtering result:
reads passed filter: 17411386
reads failed due to low quality: 6896
reads failed due to too many N: 486
reads failed due to too short: 254
reads with adapter trimmed: 17183714
bases trimmed due to adapters: 954390822

Duplication rate (may be overestimated since this is SE data): 20.678%

JSON report: fastp.json
HTML report: fastp.html

fastp -i H1_R1.fastq.gz -o H1_R1.fastp.fastq.gz 
fastp v0.23.4, time used: 66 seconds

so I tried to inspect my fastq (I highlighted my adapter with **)

@VH00643:236:AACFYWKHV:1:1101:46597:1076 1:N:0:CACTGAACCG+TCCTAAGGTT
GGTAAAGGCTCACCAAGGCGAACTGTAGGCACCATCAATTAGCCAATAGAT**AGATCGGAAGAGCACACGTCTGAACTCCAGTCA**CCACTGAACCGATCTCG
+
CC;---CCCC;CC;;CCCCCC;CCCCCCCCC;C-C-CCCC-C-CCCC-C-C;C-CCC;CCC;CCC;CCCCCCCCCCC;CCCCCC-C-CCC-;;CCC;-;-C
@VH00643:236:AACFYWKHV:1:1101:48680:1151 1:N:0:CACTGAACCG+TCCTAAGGTT
GTAGGAGTCTTGCGCAACTGTAGGCACCATCAATGCCACATAAAAA**AGATCGGAAGAGCACACGTCTGAACTCCAGTCA**CCACTGAACCGATGTGGGTTGG
+
;--CC-C--CCCCCCCCC;CCCCCCC;-CCC-CCC-C;C-C--;;CCC;CCCCCCC-CCC;-;CC-CCCCCCC--CCC;C---CC;-C-CCC--CCC-C-C
@VH00643:236:AACFYWKHV:1:1101:44192:1208 1:N:0:CACTGAACCG+TCCTAAGGTT
GTGGCTAGAGGAACTGTAGGCACCATCAATTCTTGCCAACCT**AGATCGGAAGAGCACACGTCTGAACTCCAGTCA**CCACTGAACCGATCTGGGTTGGGGTG
+
;-CC-C-C-CC;CCCCC;CC;CCC;C-CCCCCCCCCCCCC-CCC;CCCCCCC;CCCCCCCCCCCCCCC;CCCCCCCC;CCC;CC-CCC---C-CC-C-C-C
@VH00643:236:AACFYWKHV:1:1101:48661:1208 1:N:0:CACTGAACCG+TCCTAAGGTT
TGAAGCCAGAGAACTGTAGGCACCATCAATTTCCAGCGCTAC**AGATCGGAAGAGCACACGTCTGAACTCCAGTCA**CCACTGAACCGATCTGGGTTTGCGTG
+
-C--C;;-C-C-C;CCC-CCCCCC-CCCCCCCCCCCCC;CC;-CCCCCCCCC;C;-;CCCCCCC;CCCCC-CCCCCC;CCC;;C-C-C---C--C-C-C;C
@VH00643:236:AACFYWKHV:1:1101:42582:1227 1:N:0:CACTGAACCG+TCCTAAGGTT
TTCAGTGGGGTGGCGAACTGTAGGCACCATCAATCCCGATTCCCGA**AGATCGGAAGAGCACACGTCTGAACTCCAGTCA**CCACTGAACCGATCTCGGTTTG
+
-C;;CCCCCCCCCCCC-CCCCCCCCCCC;CCCCC;CCCCCCCC;CCCCCC;;CCCC-C-;;CCCCC;CCCCCCCCCC;CCC-CCCC;;CCCCCC-C--CC;
@VH00643:236:AACFYWKHV:1:1101:48869:1227 1:N:0:CACTGAACCG+TCCTAAGGTT
GGTGGAGTGATTTGTCTGGTTAATTCCGTTAACGAACGAGACCTTAACCTAACTGTAGGCACCATCAATAACATCGAGGTGTGTTCGGAAGTGCACACGTC
+
;C;CC-CCC-CCCCCCCCCCC;-CCCCCCC-C;C;CCCCC;--CCCC--C;;CCCCCCC-;-C-CC-CC--;;CCC--C-C-C-C-CCC-C-C----;CC-
@VH00643:236:AACFYWKHV:1:1101:56973:1227 1:N:0:CACTGAACCG+TCCTAAGGTT
GCTAAGCAGGGTCGGGCCTAACTGTAGGCACCATCAATGACGCATAACTA**AGATCGGAAGAGCACACGTCTGAACTCCAGTCA**CCACTGAACCGATCTGGG
+
;-C--C--CCCCCCCCCCCCCCCCC-CCCCCCCCCCC-C;CCCCCCCCCCCC-C;;;CCC-CCCCCCCCC;CCCCCCCCCCC;CCCC;CC-C;CCC;C-C-
@VH00643:236:AACFYWKHV:1:1101:45479:1246 1:N:0:CACTGAACCG+TCCTAAGGTT
GTGGGAGGGGCCGTAACTGTAGCACCATCAATAGACATTAAGCA**AGATCGGAAGAGCACACGTCTGAACTCCAGTCA**CCACTGAACGGATGTGGGTTGGGG
+
C;CCCCCCCCCCCC-CCCCCCCCCCCCCCCCC-C-CCC;CCCCCCC-CCCCCCCCCCC-C-CCCCCCCCCCCCCCC;-;-CCCCC--C-C--CC-CC-;;C
@VH00643:236:AACFYWKHV:1:1101:56841:1246 1:N:0:CACTGAACCG+TCCTAAGGTT
GGGGCTTGGAAGAATCAGCGGGGAACTGTAGGCACCATCAATAGCACAGCATGAAGATCGGAAGAGCACACGTCTGAACTCCTGTCACCACTGAACGGATG
+
-CCC-;CCC-CC-CCC-CCCCCCCCCCCCCCCC;CCCCCCC-;CCCCCCC-CCCCC-CCCCCCCCC-CCCCCCCCCCCCCCC-CC--;C-CCC;---C-C;
@VH00643:236:AACFYWKHV:1:1101:59719:1284 1:N:0:CACTGAACCG+TCCTAAGGTT
GTGTCCTCCGGGAAGTCCTGAACTGTAGGCACCATCAATGGCCTCTACTCG**AGATCGGAAGAGCACACGTCTGAACTCCAGTCA**CCACTGAACCGATGTGG
+
C-C-;;CC-CCC;CCCCCCCCCCCCC-CCC;CCCCC-CCCCCCCCCCCC;CCC-C;CCCCCCC;;CC;CCCCCCCCCC-CCCCC;C-CCCCCC-CCC---C
@VH00643:236:AACFYWKHV:1:1101:35690:1303 1:N:0:CACTGAACCG+TCCTAAGGTT

But I think there is also that needs to be filtered because the mean length of peaks is 46 bp.

[Chironex]

dear genomax, they replied "We have used the QIAseq miRNA UDI Library Kit. The used UDIs you can be found in the report we have provided you in your sFTP folder in Article II." Where I have to add this adapters? 5' end?

[GenoMax]

Looks like there may be more than one kit version. You should be able to find the manuals from Qiagen's site. Here is one version: https://genome.med.harvard.edu/documents/libraryPrep/QiagensmRNAmiRNAProtocol.pdf On page 57 you can find the adapter sequences. Generally adapters are ligated to 3'-end of miRNA's. But there may be other library specific adapters on 5'-end.

[Chironex]

Ok, thank you very much!!! I used the adapted indicated here and by your file and I was able to get a decent number of reads with a peak between 16 and 23 bp. I will go forward with other analysis and I will see if there is a decent alignment rate!