Multimapping 21-mers in rde1Δ

In the tutorial Cryptococcus siRNAs it was observed that in the strain lacking Rde1 most multimappers had a length of 21 nt, rather than 22 nt as found in other strains. Appearance of 21-mer multimappers can be analyzed further:

To find counts for reads with matches from 20 to 23 nucleotides in the alignments, as shown in Table 1, the following approach can be taken. As a control, cDNAs of multimappers are counted.

First, change directory, for the multimapping cDNAs in wild type:

• cd STAR-analysis0-h99_collapsed/SRR646636_collapsed_0mismatch-h99/

Then, get multimapper cDNA counts with samtools and grep; loop through a series of lengths:

• for i in 0 1 2 3; do samtools view samtoolsAligned.sortedByCoord.out.bam | grep -v "NH:i:1" | grep -cw "2${i}M"; done (-v for take the inverse, -c for count and -w for search whole word).

Table 1

Sample		cDNA counts (collapsed)				Multimapper read counts (uncollapsed)
Name	SRA	20-mer	21-mer	22-mer	23-mer	20-mer	21-mer	22-mer	23-mer
rde1Δ	SRR8697586	9973	12924	15972	12859	375984	723264	550248	446489
	SRR8697587	9627	13294	16094	13380	309397	593993	451772	403171
wild type	SRR646636	32990	45271	53001	42888	322217	741139	1557254	1177411
rde2Δ	SRR8697588	24192	34445	41074	34790	603920	1264609	1709966	1291535
	SRR8697589	16108	31314	43004	37049	268841	865740	1345887	1053252

Repeat the command for the mutants after:

• cd ../../STAR-analysis0-h99_collapsed/SRR8697586_collapsed_0mismatch-h99/

• cd ../../STAR-analysis0-h99_collapsed/SRR8697587_collapsed_0mismatch-h99/

and for their uncollapsed, multimapping reads:

• cd ../../STAR-analysis0-h99_uncollapsed/SRR8697586_uncollapsed_0mismatch-h99/

• cd ../../STAR-analysis0-h99_uncollapsed/SRR8697587_uncollapsed_0mismatch-h99/

or those of the wild-type, or rde2Δ as controls:

• cd ../../STAR-anaΔ lysis0-h99_uncollapsed/SRR646636_uncollapsed_0mismatch-h99/

• cd ../../STAR-analysis0-h99_collapsed/SRR8697588_collapsed_0mismatch-h99/

• cd ../../STAR-analysis0-h99_collapsed/SRR8697589_collapsed_0mismatch-h99/

• cd ../../STAR-analysis0-h99_uncollapsed/SRR8697588_uncollapsed_0mismatch-h99/

• cd ../../STAR-analysis0-h99_uncollapsed/SRR8697589_uncollapsed_0mismatch-h99/

Although the number of different multimapping reads (cDNAs) peaks for 22-mers in each strain; the counts are highest for 21-mers in rd1Δ.

To get the numbers for 21-mers for each chromosome do (in each directory with uncollapsed reads):

• for i in {1..14}; do echo "$i: $(samtools view samtoolsAligned.sortedByCoord.out.bam | grep -v "NH:i:1" | awk "\$3==$i" | grep -cw "21M" )"; done

(awk '$3=1' will select lines when the third field - the chromosome name - is 1; in the full command the chromosome names are passed in as $i, a bash-variable; because bash has to process the awk command as text, the awk-qualifier $3 needs escaping and the whole is in between double quotes to allow interpretation of $i; single quotes work when comparing to a number, instead of a variable value.)

Table 2 collects the results for uncollapsed data:

Table 2

Sample	rde1Δ		wildtype	rde2Δ
SRA	SRR8697586	SRR8697587	SRR646636	SRR8697588	SRR8697589
Chromosome	Multimapper read counts for 21-mer siRNAs
1	11178	11494	99228	55428	38513
2	687818	558685	130595	922071	635397
3	2623	2736	59920	25222	16862
4	1286	1276	28462	7730	6217
5	1089	1115	69856	26562	16624
6	891	849	63288	30566	21108
7	3737	3979	26993	24209	15618
8	784	846	28859	16604	10888
9	7773	7514	67742	52420	37142
10	526	554	15952	13336	7548
11	1491	1382	72158	36953	23800
12	1000	841	13926	8191	5484
13	1227	1299	50439	38230	25546
14	491	452	11011	4850	3442

With the exception of chromosome 2, 8 to 30-fold more reads are mapped to all chromosomes of the wild type than in rde1Δ strains. Apart from chromosome 2, also less reads are found for chromosomes in rde2Δ compared to wildtype. In the rde1Δ mutants, the vast majority of reads map to chromosome 2. Because the counted reads represent a mix of specific and unspecific reads of 21 nt, hits might be mapped to unspecific regions. Given that rRNA reads are multimappers on chromosome 2, these numbers could point to the rDNA locus on this chromosome.

With the following commands we can narrow this down; run in each folder SRR..._uncollapsed_0mismatch-h99/ a command to find the number of reads that map to chromosome 2:

• samtools view samtoolsAligned.sortedByCoord.out.bam | grep -v "NH:i:1" | awk "\$3==2" | grep -cw "21M"

and then one that retrieves the counts of reads that align within the rDNA locus between 2:270000-284000:

• samtools view samtoolsAligned.sortedByCoord.out.bam | grep -v "NH:i:1" | awk "\$3==2" | awk '$4 > 270000 && $4 < 284000' | grep -cw "21M"

for the totals you could run the command used for Table 1:

• samtools view samtoolsAligned.sortedByCoord.out.bam | grep -v "NH:i:1" | grep -cw "21M"

This gives :

Table 3

Sample	rde1Δ		wildtype	rde2Δ
SRA	SRR8697586	SRR8697587	SRR646636	SRR8697588	SRR8697589
Chromosome	Multimapper read counts for 21-mer siRNAs
Total	723264	593993	738429	1264609	865740
2:start-end	687818	558685	130595	922071	635397
2:270000-284000	684210	555345	32329	882954	611409

This approach points out that most 21-mers in rde1Δ are linked to the rDNA. Whether these represent siRNAs or rRFs is not clear. An analysis with coalispr region, however, revealed that siRNAs antisense to rRNA can be generated which were highly abundant in rde1Δ, rde2Δ (mostly against LSU rRNA) or rrp6Δ (mostly SSU rRNA; see Figure 13 in the tutorial).

We gain insight of the number of different multimappers by counting cDNAs with 21N gaps, and by checking the mapping strand (using the Flag), where the derive from:

• cd ../../STAR-analysis0-h99_collapsed/

Run similar commands as above for the uncollapsed data and get:

• total number of cDNAs for multimapping 21-mers:

  samtools view samtoolsAligned.sortedByCoord.out.bam | grep -v "NH:i:1" | grep -cw "21M"

• total number of cDNAs for multimapping 21-mers on chromosome 2:

samtools view samtoolsAligned.sortedByCoord.out.bam | grep -v "NH:i:1" | awk '$3==2' | grep -cw "21M"

• total number of cDNAs for multimapping 21-mers at the 2 rDNA loci on chromosome 2:

  samtools view samtoolsAligned.sortedByCoord.out.bam | grep -v "NH:i:1" | awk '$3==2' | awk '$4 > 270000 && $4 < 284000' | grep -csw "21M"

• of those on the positive strand (+), FLAG 0:

  samtools view samtoolsAligned.sortedByCoord.out.bam | grep -v "NH:i:1" | awk '$3==2' | awk '$2==0' | awk '$4 > 270000 && $4 < 284000' | grep -csw "21M"

• or negative strand (-), from (pre-)rRNA, FLAG 16:

  samtools view samtoolsAligned.sortedByCoord.out.bam | grep -v "NH:i:1" | awk '$3==2' | awk '$2==16' | awk '$4 > 270000 && $4 < 284000' | grep -csw "21M"

• or as secondary alignment on the positive strand (++), FLAG 256:

  samtools view samtoolsAligned.sortedByCoord.out.bam | grep -v "NH:i:1" | awk '$3==2' | awk '$2==256' | awk '$4 > 270000 && $4 < 284000' | grep -csw "21M"

• or as secondary alignment on the negative strand (--), from (pre-)rRNA, FLAG 272:

  samtools view samtoolsAligned.sortedByCoord.out.bam | grep -v "NH:i:1" | awk '$3==2' | awk '$2==272' | awk '$4 > 270000 && $4 < 284000' | grep -csw "21M"

The outcome, shown below in bold, shows that, compared to wild type, a relatively large number of different 21-mer siRNAs are generated in the absence of Rde1 or Rde2 against rRNA sequences or rRFs.

Table 4

Sample	rde1Δ		wildtype	rde2Δ
SRA	SRR8697586	SRR8697587	SRR646636	SRR8697588	SRR8697589
Chromosome	Multimapper cDNA counts for 21-mer reads
Total	12406	13294	45271	34445	31314
2:start-end	4871	4730	7187	7891	4045
2:270000-284000	3984	3609	2328	4425	7247
+  FLAG 0	668	514	35	856	762
-  FLAG 16	3228	3022	2231	3472	3203
++ FLAG 256	4	1	4	8	2
-- FLAG 272	84	72	58	89	78

After these tests return to Burke-2019/:

• cd ../../