1. Staphylococcus aureus, Vibrio anguillarum and Vibrio alginolyticus (these are Branchiostoma pathogens according to Huang et al., 2014 and Zou et al., 2016) genomes were downloaded from the NCBI "Genome" database on November 12, 2018. Genome sequence files are named 'Staphylococcus_aureus_genome.fa.bz2', 'Vibrio_anguillarum_genome.fa.bz2' and 'Vibrio_alginolyticus_genome.fa.bz2'. Mapping non-Branchiostoma reads on these genomes:

``for species in Staphylococcus_aureus Vibrio_alginolyticus Vibrio_anguillarum;do bowtie2-build $species'_genome.fa' $species'_genome';done;for sample in embryon_15h embryon_36h embryon_60h embryon_8h femelle male;do for lib in `seq 1 4`;do for species in Staphylococcus_aureus Vibrio_alginolyticus Vibrio_anguillarum;do bowtie2 -x $species'_genome' -U ../Transcriptome/Extragenomic_extratranscriptomic_reads_$sample'_'$lib'.fastq' --no-unal --quiet -S $species'_mapping_'$sample'_'$lib'.sam';done;done;done;for species in Staphylococcus_aureus Vibrio_alginolyticus Vibrio_anguillarum;do for sample in embryon_15h embryon_36h embryon_60h embryon_8h femelle male; do for lib in `seq 1 4`; do echo "Orientation Position Sequence" > $species'_genome_mappers_in_'$sample'_'$lib'.dat';sed '1,/^@PG\t/ d' $species'_mapping_'$sample'_'$lib'.sam' | awk '{print $2,$4,$10}' >> $species'_genome_mappers_in_'$sample'_'$lib'.dat';done;done;done;for species in Staphylococcus_aureus Vibrio_alginolyticus Vibrio_anguillarum;do genome_size=`~/Fuses_lines_clean.pl $species'_genome.fa' | tail -1 | wc -m`;genome_size=`echo $genome_size"-1" | bc`;Rscript R_commands_mapping_on_bacterial_genomes $species $genome_size;done``

Counting the total number of reads for each of these 3 genomes (when focusing on reads that do not match either the Branchiostoma genome or the Branchiostoma transcriptome):

``cat Extragenomic_extratranscriptomic_reads_matching_Staphylococcus_aureus_in_* | grep -c '>' # Result: 1581221``
``cat Extragenomic_extratranscriptomic_reads_matching_Vibrio_alginolyticus_in_* | grep -c '>' # Result: 810894``
``cat Extragenomic_extratranscriptomic_reads_matching_Vibrio_anguillarum_in_* | grep -c '>' # Result: 780138``
``cat list_S_aureus list_V_alginolyticus list_V_alginolyticus list_V_anguillarum list_V_anguillarum list_V_anguillarum list_V_anguillarum | sort | uniq -c | awk '{print $1}' | sort | uniq -c # Result: 1457122 unique to S. aureus; 113398 unique to V. alginolyticus; 25048 shared only by S. aureus and V. alginolyticus; 103153 unique to V. anguillarum; 4537 shared only by S. aureus and V. anguillarum; 577934 ahread only by V. alginolyticus and V. anguillarum; 94514 shared by all three species``

Size distribution histograms and logos for the small RNAs that map on these pathogenic genomes (Supplementary File S1 sections 7, 8 and 9):

``sample=embryon_8h;nohup ./Chef.sh $sample > nohup_$sample'.out'&;sample=embryon_15h;nohup ./Chef.sh $sample > nohup_$sample'.out'&;sample=embryon_36h;nohup ./Chef.sh $sample > nohup_$sample'.out'&;sample=embryon_60h;nohup ./Chef.sh $sample > nohup_$sample'.out'&;sample=femelle;nohup ./Chef.sh $sample > nohup_$sample'.out'&;sample=male;nohup ./Chef.sh $sample > nohup_$sample'.out'&``


2. Assembling extragenomic, extranscriptomic reads into "foreign DNA" contigs, using velvet (with k varying from 9 to 19 as in pipeline "S" described in Massart et al., 2018, and informing Velvet that sequencing is strand-specific):

``cat Extragenomic_extratranscriptomic_reads_*.fastq > All_libraries.fastq;velveth All_libraries_k_9_19 9,21,2 -fastq -strand_specific All_libraries.fastq;for k in `seq 9 2 19`;do velvetg All_libraries_k_9_19_$k -min_contig_lgth 26;done``

Selecting contigs at least 50 bp long, and preparing several subsets based on their lengths (50-54 bp long; 55-59 bp long; ...; 95-99 bp long; 100 bp or longer):

``for k in `seq 9 2 19`;do ~/Fuses_lines_clean.pl All_libraries_k_9_19_$k/contigs.fa > Fused_contigs_k_$k'.fa';done;cat Fused_contigs_k_* | grep -B 1 '^[ACGT]\{100\}' | grep -v '\-\-' > long_contigs.fa;for i in `seq 0 9`;do start=`echo $i"*5+50" | bc`;end=`echo $i"*5+54" | bc`;cat Fused_contigs_k_* | grep -B 1 '^[ACGT]\{'$start,$end'\}$'  | grep -v '\-\-' > medium_contigs_$i'.fa';done``

Contig annotation: BLAST'ing these 'long_contigs.fa' and 'medium_contigs_*.fa' files on the NCBI web server on October 31, 2018; database: "Nucleotide collection (nr/nt)", Program  Selection: "Optimize for Highly similar sequences (megablast)". Results downloaded as 11 XML files ('XK9RHVPE01R-Alignment_Single_file.xml' for contigs 100 bp or longer; 'medium_0-Alignment.xml' to 'medium_9-Alignment.xml' for 50-99 bp-long contigs). Then:

``cat medium_*.xml > All_medium_Alignment.xml;cat XK9RHVPE01R-Alignment_Single_file.xml All_medium_Alignment.xml > All_medium_and_long_contigs.xml;./Module_parse_XML.pl All_medium_and_long_contigs.xml;awk -F '\t' '{print $2}' Annotated_best_hits_from_All_medium_and_long_contigs.txt | sort | uniq -c | sort -g``

(of note: these contigs are not unified: some of them were captured mutliply, with various values for the "k" parameter in the velvet analysis)


3. For the preparation of Table 2: classification of phylogenetically-annotated contigs:

``grep -v Branchiostoma Annotated_best_hits_from_All_medium_and_long_contigs.txt | awk -F '\t' '{print $2}' | sed 's|/|\
|g' | sort | uniq > rough_species_list``

(this is a "rough" list, because it assumes that the first 2 words of sequence annotation are the genus + species names; that's most often the case, but some manual correction will be needed afterwards for exceptions)

``n=`grep -vc '^[0-9]' rough_species_list`;for i in `seq 1 $n`;do genus=`grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | awk '{print $1}'`;species=`grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | awk '{print $2}'`;echo "genus="$genus" species="$species;sleep 0.2;wget -O esearch_$i http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=taxonomy\&term=$genus'+'$species;id=`grep '^<Id>' esearch_$i | sed -e 's|^<Id>||' -e 's|</Id>||'`;sleep 0.2;wget -O efetch_$i http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=taxonomy\&id=$id;grep '^ *<Lineage>' efetch_$i | sed -e 's|^ *<Lineage>||' -e 's|</Lineage>$||' > Phylogenetic_annotation_$i;done``

845 values of $i have an empty 'Phylogenetic_annotation_$i' file; for 198 of them, that's because the species name was not given (e.g., names like "Anaeromyxobacter sp.", with "sp." instead of the species name). So for these ones, let's just use the genus name:

``for i in `wc -l Phylogenetic_annotation_* | grep '^ *0 ' | sed 's|.*_||'`;do genus=`grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | awk '{print $1}'`;species=`grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | awk '{print $2}'`;if test "$species" = "sp.";then sleep 0.2;wget -O esearch_$i http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=taxonomy\&term=$genus;id=`grep '^<Id>' esearch_$i | sed -e 's|^<Id>||' -e 's|</Id>||'`;sleep 0.2;wget -O efetch_$i http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=taxonomy\&id=$id;grep '^ *<Lineage>' efetch_$i | sed -e 's|^ *<Lineage>||' -e 's|</Lineage>$||' > Phylogenetic_annotation_$i;fi;done``

For 19 of the remaining 647 empty phylogenies, the genus name was given between square brackets. Correcting that:

``for i in `wc -l Phylogenetic_annotation_* | grep '^ *0 ' | sed 's|.*_||'`;do genus=`grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | awk '{print $1}'`;species=`grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | awk '{print $2}'`;if test `echo $genus | grep -c '^\['` -eq 1;then sleep 0.2;wget -O esearch_$i http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=taxonomy\&term=$genus;id=`grep '^<Id>' esearch_$i | sed -e 's|^<Id>||' -e 's|</Id>||'`;sleep 0.2;wget -O efetch_$i http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=taxonomy\&id=$id;grep '^ *<Lineage>' efetch_$i | sed -e 's|^ *<Lineage>||' -e 's|</Lineage>$||' > Phylogenetic_annotation_$i;fi;done``

For 128 of the remaining 628 empty phylogenies, the genus name was preceded by "PREDICTED: ". Correcting that:

``for i in `wc -l Phylogenetic_annotation_* | grep '^ *0 ' | sed 's|.*_||'`;do genus=`grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | awk '{print $1}'`;species=`grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | awk '{print $2}'`;if test `echo $genus | grep -c '^PREDICTED:'` -eq 1;then sleep 0.2;wget -O esearch_$i http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=taxonomy\&term=$species;id=`grep '^<Id>' esearch_$i | sed -e 's|^<Id>||' -e 's|</Id>||'`;sleep 0.2;wget -O efetch_$i http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=taxonomy\&id=$id;grep '^ *<Lineage>' efetch_$i | sed -e 's|^ *<Lineage>||' -e 's|</Lineage>$||' > Phylogenetic_annotation_$i;fi;done``

For 88 of the remaining 500 empty phylogenies, there is a dot in the (normally) genus name (in many cases that's because the genus was abbreviated, and concatenated without space with the species name, e.g., "D.melanogaster"). Correcting that (N.B.: for some of them, the species name was too ambiguous for a clear attribution: we will not correct those):

``for i in `wc -l Phylogenetic_annotation_* | grep '^ *0 ' | sed 's|.*_||'`;do species="";genus=`grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | awk '{print $1}'`;case "$genus" in "S.cerevisiae") genus=Saccharomyces;species=cerevisiae;;"D.melanogaster") genus=Drosophila;species=melanogaster;;"L.lactis") genus=Lactococcus;species=lactis;;"H.sapiens") genus=Homo;species=sapiens;;"D.yakuba") genus=Drosophila;species=yakuba;;"M.musculus") genus=Mus;species=musculus;;"d.melanogaster") genus=Drosophila;species=melanogaster;;"D.erecta") genus=Drosophila;species=erecta;;"B.floridae") genus=Branchiostoma;species=floridae;;"A.pisum") genus=Acyrthosiphon;species=pisum;;esac;if test "$species" != "";then sleep 0.2;wget -O esearch_$i http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=taxonomy\&term=$genus'+'$species;id=`grep '^<Id>' esearch_$i | sed -e 's|^<Id>||' -e 's|</Id>||'`;sleep 0.2;wget -O efetch_$i http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=taxonomy\&id=$id;grep '^ *<Lineage>' efetch_$i | sed -e 's|^ *<Lineage>||' -e 's|</Lineage>$||' > Phylogenetic_annotation_$i;fi;done``

For some of the remaining 461 empty phylogenies, some keywords ("cerevisiae", "Drosophila", ...) appeared in either the (normally) "genus" or "species". Correcting for that (only for keywords appearing at least twice in the list: the other ones will be considered as "poorly annotated" and won't be corrected manually):

``for i in `wc -l Phylogenetic_annotation_* | grep '^ *0 ' | sed 's|.*_||'`;do clade='';line=`grep -v '^[0-9]' rough_species_list | head -$i | tail -1`;if test `grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | grep -c 'S.cerevisiae'` -ne 0;then clade=Saccharomyces;fi;if test `grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | grep -c '[Yy]east'` -ne 0;then clade=Ascomycota;fi;if test `grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | grep -c 'Bacterium'` -ne 0;then clade=Bacteria;fi;if test `grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | grep -c 'Drosophila'` -ne 0;then clade=Drosophila;fi;if test `grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | grep -c 'Human'` -ne 0;then clade=Homo;fi;if test `grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | grep -c 'L.lactis'` -ne 0;then clade=Lactococcus;fi;if test `grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | grep -c 'H.sapiens'` -ne 0;then clade=Homo;fi;if test `grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | grep -c 'Bacillus'` -ne 0;then clade=Bacilli;fi;if test `grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | grep -c 'Vibrio'` -ne 0;then clade=Vibrio;fi;if test `grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | grep -c 'Rickettsia'` -ne 0;then clade=Rickettsia;fi;if test `grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | grep -c 'Rat'` -ne 0;then clade=Rattus;fi;if test `grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | grep -c 'Galactomyces'` -ne 0;then clade=Galactomyces;fi;if test `grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | grep -c 'Escherichia'` -ne 0;then clade=Escherichia;fi;if test `grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | grep -c 'D.yakuba'` -ne 0;then clade=Drosophila;fi;if test `grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | grep -c 'Acinetobacter'` -ne 0;then clade=Acinetobacter;fi;if test `grep -v '^[0-9]' rough_species_list | head -$i | tail -1 | grep -c 'Leuconostoc'` -ne 0;then clade=Leuconostoc;fi;if test "$clade" != "";then sleep 0.2;wget -O esearch_$i http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=taxonomy\&term=$clade;id=`grep '^<Id>' esearch_$i | sed -e 's|^<Id>||' -e 's|</Id>||'`;sleep 0.2;wget -O efetch_$i http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=taxonomy\&id=$id;grep '^ *<Lineage>' efetch_$i | sed -e 's|^ *<Lineage>||' -e 's|</Lineage>$||' > Phylogenetic_annotation_$i;fi;done``

Using that manually corrected list, annotating BLAST hits for a few clades of interest:

``n=`cat rough_species_list | wc -l`;j=0;for i in `seq 1 $n`;do if test `head -$i rough_species_list | tail -1 | grep -c '^[0-9]'` -eq 0;then j=`echo $j"+1" | bc`;head -$i rough_species_list | tail -1 | sed 's|$| phylogeny_'$j'|';else head -$i rough_species_list | tail -1 | sed 's|$| phylogeny_NA|';fi;done > Phylogenic_annotation_of_rough_species_list;clades_of_interest='Bacteria Fungi Viridiplantae Metazoa Branchiostoma Viruses';n=`cat Annotated_best_hits_from_All_medium_and_long_contigs.txt | wc -l`;head -1 Annotated_best_hits_from_All_medium_and_long_contigs.txt | sed 's|$|\tClade|' > Phylogenetically_Annotated_best_hits_from_All_medium_and_long_contigs.txt;for i in `seq 2 $n`;do phyl='';annot=`head -$i Annotated_best_hits_from_All_medium_and_long_contigs.txt | tail -1 | awk -F '\t' '{print $2}'`;echo $annot | sed 's|/|\
|g' | awk '{print $1,$2}' > tmp_split;m=`cat tmp_split | wc -l`;for j in `seq 1 $m`;do term1=`head -$j tmp_split | tail -1 | awk '{print $1}' | sed -e 's|\[||g' -e 's|\]||g'`;term2=`head -$j tmp_split | tail -1 | awk '{print $2}' | sed -e 's|\[||g' -e 's|\]||g'`;index=`grep '^'$term1' '$term2 Phylogenic_annotation_of_rough_species_list | sed 's|.*phylogeny_||'`;if test "$term1" = "Branchiostoma";then phyl=`echo $phyl"/Branchiostoma"`;else for ind in `echo $index`;do if test "$ind" != "NA";then for clade in `echo $clades_of_interest`;do if test `grep -wc $clade Phylogenetic_annotation_$ind` -ne 0;then phyl=`echo $phyl"/"$clade`;fi;done;fi;done;fi;done;phyl=`echo $phyl | sed 's|^/||'`;if test "$phyl" = "";then phyl=NA;fi;head -$i Annotated_best_hits_from_All_medium_and_long_contigs.txt | tail -1 | sed 's|$|\t'$phyl'|';done >> Phylogenetically_Annotated_best_hits_from_All_medium_and_long_contigs.txt;cp Phylogenetically_Annotated_best_hits_from_All_medium_and_long_contigs.txt tmp_annot;apres=`md5sum tmp_annot | awk '{print $1}'`;avant=0;i=1;while test "$avant" != "$apres";do avant=`md5sum tmp_annot | awk '{print $1}'`;echo "iteration "$i;i=`echo $i"+1" | bc`;sed -e 's|Bacteria/Bacteria|Bacteria|g' -e 's|Fungi/Fungi|Fungi|g' -e 's|Viridiplantae/Viridiplantae|Viridiplantae|g' -e 's|Metazoa/Metazoa|Metazoa|g' -e 's|Branchiostoma/Branchiostoma|Branchiostoma|g' -e 's|Viruses/Viruses|Viruses|g' tmp_annot > tmp2_annot;mv -f tmp2_annot tmp_annot;apres=`md5sum tmp_annot | awk '{print $1}'`;done;mv -f tmp_annot Phylogenetically_Annotated_best_hits_from_All_medium_and_long_contigs.txt``

Problem: some viruses are annotated as "Blah blah virus", with "Blah blah" being the name of a species (so my pipeline annotates these contigs with the name of that host species). It has to be corrected by hand (N.B.: not for "Guanarito mammarenavirus", which is suspiciously similar to metazoan rRNA: let's rather stick to the "Metazoan" annotation for that contig):

``for contig in `grep -i virus Detailed_annotated_best_hits_from_All_medium_and_long_contigs.txt  | grep -v 'Guanarito mammarenavirus' | grep 'Escherichia virus' | awk '{print $1}'` NODE_1627332_length_36_cov_1.222222 NODE_923783_length_43_cov_7.953488 NODE_879336_length_52_cov_1.134615;do sed -i '/^'$contig'\t/ s|\(.*\)\t.*|\1\tProkaryotic_virus|' Phylogenetically_Annotated_best_hits_from_All_medium_and_long_contigs.txt;done;for contig in NODE_926073_length_49_cov_3.204082 NODE_926566_length_40_cov_2.700000 NODE_704291_length_38_cov_2.552632 NODE_980537_length_50_cov_3.280000 NODE_722627_length_48_cov_3.000000 NODE_685849_length_52_cov_11.442307;do sed -i '/^'$contig'\t/ s|\(.*\)\t.*|\1\tEukaryotic_virus|' Phylogenetically_Annotated_best_hits_from_All_medium_and_long_contigs.txt;done``

Statistics on the total contig set:

``cat medium_contigs_*.fa long_contigs.fa | grep -v '>' | sort | uniq | wc -l;for seq in `cat medium_contigs_*.fa long_contigs.fa | grep -v '>'`;do echo -n $seq | wc -m;done > tmp.dat;R CMD BATCH R_commands_contig_length_statistics;echo `head -2 out.csv | tail -1 | awk -F ',' '{print $2}'`" ($\pm$"`tail -1 out.csv | awk -F ',' '{print $2}'`")"``

Statistics on various subsets:

``for clade in Bacteria Fungi Viridiplantae Metazoa Branchiostoma Prokaryotic_virus Eukaryotic_virus;do echo $clade":";for contig in `awk -F '\t' '$4=="'$clade'" {print $1}' Phylogenetically_Annotated_best_hits_from_All_medium_and_long_contigs.txt`;do cat medium_contigs_*.fa long_contigs.fa | grep -A 1 '^>'$contig'$' | grep -v '>';done | sort | uniq > tmp.txt;for seq in `cat tmp.txt`;do echo -n $seq | wc -m;done > tmp.dat;R CMD BATCH R_commands_contig_length_statistics;echo `cat tmp.dat | wc -l`" unique contigs";echo `head -2 out.csv | tail -1 | awk -F ',' '{print $2}'`" ($\pm$"`tail -1 out.csv | awk -F ',' '{print $2}'`")";echo "";done``

Ambiguous, poorly described and others:

``cat medium_contigs_*.fa long_contigs.fa | grep -v '>' | sort | uniq > list_1;for clade in Bacteria Fungi Viridiplantae Metazoa Branchiostoma Prokaryotic_virus Eukaryotic_virus;do for contig in `awk -F '\t' '$4=="'$clade'" {print $1}' Phylogenetically_Annotated_best_hits_from_All_medium_and_long_contigs.txt`;do cat medium_contigs_*.fa long_contigs.fa | grep -A 1 '^>'$contig'$' | grep -v '>';done;done | sort | uniq > list_2;for seq in `cat list_1 list_2 list_2 | sort | uniq -c | grep '^ *1 ' | awk '{print $2}'`;do echo -n $seq | wc -m;done > tmp.dat;R CMD BATCH R_commands_contig_length_statistics;echo `cat tmp.dat | wc -l`" unique contigs";echo `head -2 out.csv | tail -1 | awk -F ',' '{print $2}'`" ($\pm$"`tail -1 out.csv | awk -F ',' '{print $2}'`")"``


4. Identification of eukaryotic virus-derived contigs:

``for contig in `sort Detailed_annotated_best_hits_from_All_medium_and_long_contigs.txt | uniq -c | sort -g | grep -v 'Escherichia virus' | grep -v ribosom | grep -i virus | awk '{print $2}'`;do cat Fused_contigs_k_*.fa | grep -A 1 '^>'$contig'$';done > Viral_contigs.fa``

(of note: these contigs are not unified: some of them were captured mutliply, with various values for the "k" parameter in the velvet analysis; out of the 6 sequences in 'Viral_contigs.fa', there are only 4 unique sequences)

``bowtie2-build Viral_contigs.fa Viral_contigs;for sample in embryon_15h embryon_36h embryon_60h embryon_8h femelle male;do for lib in `seq 1 4`;do bowtie2 -x Viral_contigs -U ../Transcriptome/Extragenomic_extratranscriptomic_reads_$sample'_'$lib'.fastq' --no-unal --quiet -S Viral_contig_mapping_$sample'_'$lib'.sam';done;done;for sample in embryon_15h embryon_36h embryon_60h embryon_8h femelle male; do for lib in `seq 1 4`; do sed '1,/^@PG\t/ d' Viral_contig_mapping_$sample'_'$lib'.sam' | awk '{print $3}'; done; done | sort | uniq -c | sort -g``

Result: of these 6 non-unified contigs (4 unified contigs), 5 (3 after unification) are hardly matched by anything (no more than 8 raw reads when summing all 24 libraries) but the last contig is heavily-matched (1012 raw reads when summing all 24 libraries). That's the "NODE_685849_length_52_cov_11.442307" contig, which matches the ATCV1 genome. Analysis of read length distribution, sequence bias analysis, and read coverage for reads matching that contig:

``echo "Orientation Position Sequence" > Champion_contig_mappers.dat;for sample in embryon_15h embryon_36h embryon_60h embryon_8h femelle male; do for lib in `seq 1 4`; do sed '1,/^@PG\t/ d' Viral_contig_mapping_$sample'_'$lib'.sam' | awk '$3=="NODE_685849_length_52_cov_11.442307" {print $2,$4,$10}';done;done >> Champion_contig_mappers.dat;for sample in embryon_15h embryon_36h embryon_60h embryon_8h femelle male; do for lib in `seq 1 4`;do echo "Orientation Position Sequence" > Champion_contig_mappers_in_$sample'_'$lib'.dat';sed '1,/^@PG\t/ d' Viral_contig_mapping_$sample'_'$lib'.sam' | awk '$3=="NODE_685849_length_52_cov_11.442307" {print $2,$4,$10}' >> Champion_contig_mappers_in_$sample'_'$lib'.dat';done;done;for sample in embryon_15h embryon_36h embryon_60h embryon_8h femelle male; do for lib in `seq 1 4`;do Rscript R_commands_mapping_on_viral_contig $sample $lib 70;done;done``


5. Matching reads on the whole ATCV1 genome (not just the constructed contig) (after downloading the ATCV1 genome from the NCBI "Genome" database on November 13, 2018):

``bowtie2-build ATCV1_genome.fa ATCV1_genome;for sample in embryon_15h embryon_36h embryon_60h embryon_8h femelle male;do for lib in `seq 1 4`;do bowtie2 -x ATCV1_genome -U ../Transcriptome/Extragenomic_extratranscriptomic_reads_$sample'_'$lib'.fastq' --no-unal --quiet -S ATCV1_mapping_$sample'_'$lib'.sam';done;done;for sample in embryon_15h embryon_36h embryon_60h embryon_8h femelle male;do for lib in `seq 1 4`;do echo "Orientation Position Sequence" > ATCV1_mappers_in_$sample'_'$lib'.dat';sed '1,/^@PG\t/ d' ATCV1_mapping_$sample'_'$lib'.sam' | awk '{print $2,$4,$10}' >> ATCV1_mappers_in_$sample'_'$lib'.dat';done;done;genome_length=`~/Fuses_lines_clean.pl ATCV1_genome.fa | tail -1 | wc -m`;genome_length=`echo $genome_length"-1" | bc`;Rscript R_commands_mapping_on_viral_genome ATCV1 $genome_length``

Size distribution and logo analysis for ATCV1-matching reads:

``sample=embryon_8h;nohup ./Chef_ATCV1.sh $sample > nohup_$sample'.out'&;sample=embryon_15h;nohup ./Chef_ATCV1.sh $sample > nohup_$sample'.out'&;sample=embryon_36h;nohup ./Chef_ATCV1.sh $sample > nohup_$sample'.out'&;sample=embryon_60h;nohup ./Chef_ATCV1.sh $sample > nohup_$sample'.out'&;sample=femelle;nohup ./Chef_ATCV1.sh $sample > nohup_$sample'.out'&;sample=male;nohup ./Chef_ATCV1.sh $sample > nohup_$sample'.out'&``


6. Compiling every contig (50 bp or longer), in order to map small RNAs on each of them:

``cat medium_contigs_*.fa long_contigs.fa | grep -v '^>' | sort | uniq | nl | sed 's|^ *\([0-9]*\)\t\(.*\)|>contig_\1\
\2|' > Unified_contigs.fa;for seq in `grep -v '^>' Unified_contigs.fa`;do echo $seq `cat medium_contigs_*.fa long_contigs.fa | grep -B 1 '^'$seq'$' | grep '>' | sed 's|^>||'`;done > Contig_name_correspondance.txt;bowtie2-build Unified_contigs.fa Unified_contigs;for sample in embryon_15h embryon_36h embryon_60h embryon_8h femelle male;do for lib in `seq 1 4`;do bowtie2 -x Unified_contigs -U ../Transcriptome/Extragenomic_extratranscriptomic_reads_$sample'_'$lib'.fastq' --no-unal --quiet -S Unified_contig_mapping_$sample'_'$lib'.sam';done;done``

For the most highly-matched contigs, drawing a map with read coverage (selecting contigs matched by a high number of reads: the cutoff is 200 ppm):

``CUTOFF=200;total_depth=`awk '{s+=$5-$6} END {print s}' ../Statistics.dat`;for sample in embryon_15h embryon_36h embryon_60h embryon_8h femelle male; do for lib in `seq 1 4`; do sed '1,/^@PG\t/ d' Unified_contig_mapping_$sample'_'$lib'.sam' | awk '{print $3}'; done; done | sort | uniq -c | sort -g > Highly_mapped_contigs.txt;for contig in `awk '$1>='$CUTOFF'*'$total_depth'/1000000 {print $2}' Highly_mapped_contigs.txt`;do seq=`grep -A 1 '^>'$contig'$' Unified_contigs.fa | tail -1`;list=`grep '^'$seq' ' Contig_name_correspondance.txt | sed 's|^'$seq' ||'`;display=$contig;for id in `echo $list`;do add=`grep -P '^'$id'\t' Phylogenetically_Annotated_best_hits_from_All_medium_and_long_contigs.txt | awk -F '\t' '{print $2" ("$4")"}'`;display=`echo $display $add`;done;echo $display;done > Annotated_Highly_mapped_contigs.txt;for contig in `awk '{print $1}' Annotated_Highly_mapped_contigs.txt`;do for sample in embryon_15h embryon_36h embryon_60h embryon_8h femelle male;do for lib in `seq 1 4`;do echo "Orientation Position Sequence" > Contig_$contig'_mappers_in_'$sample'_'$lib'.dat';sed '1,/^@PG\t/ d' Unified_contig_mapping_$sample'_'$lib'.sam' | awk '$3=="'$contig'" {print $2,$4,$10}' >> Contig_$contig'_mappers_in_'$sample'_'$lib'.dat';done;done;done;for contig in `awk '{print $1}' Annotated_Highly_mapped_contigs.txt`;do seq=`grep -A 1 '^>'$contig'$' Unified_contigs.fa | tail -1`;length=`echo -n $seq | wc -m`;Rscript R_commands_mapping_on_contig $contig $length;if test "$contig" = "contig_27834";then Rscript R_commands_legended_mapping_on_contig $contig $length;fi;done``

Size distribution and logo analysis for contig-matching reads, for the 5 constructed contigs with at least 200 ppm reads over the pooled 24 libraries:

``CUTOFF=200;total_depth=`awk '{s+=$5-$6} END {print s}' ../Statistics.dat`;for contig in `awk '$1>='$CUTOFF'*'$total_depth'/1000000 {print $2}' Highly_mapped_contigs.txt`;do for sample in embryon_8h embryon_15h embryon_36h embryon_60h femelle male;do ./Chef2.sh $sample $contig;sleep 20;done;done``
