Objectif

Est il possible de déterminer le sexe du veau à venir par une méthode plus facile et plus tôt dans la période de gestation.

Idée : séquençage de l'ADN circulant après une simple prise de sang et détection de séquences originaires du chromosome Y

Partenaire

Dominique Rocha, GABI

Mat & Meth

data

séquençage de 31 échantillons de vache en gestation ayant donné naissance à des mâles ou à des femelles.

séquençage 2x151pb mais avec une forte variabilité en nombre de séquence et en taille:

- entre 26 589 séquences et 840 178 séquences (moyenne à 342 051)
- entre 37 et 151pb (taille largement majoritaire)

rapport fastqc R1 : report_R1.html
rapport fastqc R2 : report_R2.html
version :FastQC_v0.11.7 MultiQC-v1.7

génome de référence à utiliser : https://www.ncbi.nlm.nih.gov/assembly/GCA_000003205.6/

méthodes

1. masquage du génome de référence

téléchargement des régions répétées : https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/003/205/GCA_000003205.6_Btau_5.0.1/GCA_000003205.6_Btau_5.0.1_rm.out.gz

module load bioinfo/bedtools2-2.29.0

zcat ref_bwa1/GCA_000003205.6_Btau_5.0.1_rm.out.gz |tail -n +4 | awk '{print $5"\t"$6"\t"$7}' > ref_bwa1/repeat.bed
bedtools maskfasta -soft -fi ref_bwa1/GCA_000003205.6_Btau_5.0.1_genomic.fna -bed ref_bwa1/repeat.bed -fo ref_bwa1/GCA_000003205.6_Btau_5.0.1_genomic.mask.fna

2. alignement

alignement des séquences sur le génome de référence masqué

module load bioinfo/bwa-0.7.17
module load bioinfo/samtools-1.10

bwa index ref_bwa1/GCA_000003205.6_Btau_5.0.1_genomic.mask.fna

mkdir bwa1

c=1
for R1 in data/*_R1_*.fastq.gz
do
sample=`echo $R1 | cut -f 1 -d "_" | cut -f 2 -d '/'`
R2=`echo $R1 | sed 's/_R1_/_R2_/'`
job_ID=`sbatch -J bwa1_$sample -o LOG/%x.out -e LOG/%x.err --cpus-per-task=8 --mem=10G --wrap="bwa mem -t8 -R '@RG\tID:$c\tSM:$sample\tPL:Illumina\tPU:PU\tLB:LB' -o bwa1/$sample.sam ref_bwa1/GCA_000003205.6_Btau_5.0.1_genomic.mask.fna $R1 $R2" | awk '{print $NF}'`
sbatch -J bwa1_seff_$sample -o LOG/%x.out -e LOG/%x.err  --dependency=afterok:$job_ID  --wrap="seff $job_ID > bwa1/$sample.seff" 
job_ID2=`sbatch -J sort_bam_$sample -o LOG/%x.out -e LOG/%x.err --cpus-per-task=8 --mem=10G --dependency=afterok:$job_ID --wrap="samtools view -Sh bwa1/$sample.sam -b | samtools sort -@ 8 -o bwa1/$sample.bam "|awk '{print $NF}'`
sbatch -J bwa1_bamIdx_$sample -o LOG/%x.out -e LOG/%x.err --dependency=afterok:$job_ID2 --wrap="samtools index bwa1/$sample.bam " 
let c=$c+1
done

flagstat pour compter le nombre de lectures alignées

for i in bwa1/*.bam
do
sample=`basename $i | sed 's/\.bam//'`
sbatch -J bwa1_flagstat_$sample -o LOG/%x.out -e LOG/%x.err --wrap="samtools flagstat $i > bwa1/`basename $i | sed 's/\.bam/\.flagstat/'`" 
sbatch -J bwa1_flagstat_chrY_$sample -o LOG/%x.out -e LOG/%x.err --wrap="samtools view -hb $i CM001061.2 | samtools flagstat - > bwa1/`basename $i | sed 's/\.bam/\.chrY_flagstat/'`" 
sbatch -J bwa1_flagstat_properly_only_chrY_$sample -o LOG/%x.out -e LOG/%x.err --wrap="samtools view -h -f2 $i CM001061.2 | grep -v XA:Z |samtools view -Shb - | samtools flagstat - > bwa1/`basename $i | sed 's/\.bam/\.properly_onlyChrY_flagstat/'`" 
sbatch -J bwa1_flagstat_properly_q30_only_chrY_$sample -o LOG/%x.out -e LOG/%x.err --wrap="samtools view -h -f2 -q 30 $i CM001061.2 | grep -v XA:Z |samtools view -Shb - | samtools flagstat - > bwa1/`basename $i | sed 's/\.bam/\.q30_properly_onlyChrY_flagstat/'`" 
done

# fichier bilan

echo -e "#sample_name\tCPU_time\tMEM_used\tnb_reads\tnb_read_mapped\tnb_read_properly_mapped\tnb_read_mapped_on_chrY\tpercent_read_mapped_on_chrY\tnb_read_properly_mapped_only_on_chrY\tpercent_read_properly_mapped_only_on_chrY\tnb_read_mapq30_properly_mapped_only_on_chrY" > bwa1/nb_aln.tsv
for i in bwa1/*.bam
do
sample=`basename $i | sed 's/\.bam//'`
cpu=`grep "CPU Utilized:" bwa1/$sample.seff | awk '{print $NF}'`
mem=`grep "Memory Utilized:" bwa1/$sample.seff | awk '{print $3" "$4}'`
tot_read=`grep read1 bwa1/$sample.flagstat | awk '{print $1*2}'`
tot_mapped=`grep -m1 mapped bwa1/$sample.flagstat | awk '{print $1}'`
tot_properly_paired=`grep "properly paired" bwa1/$sample.flagstat | awk '{print $1}'`
tot_mapped_chrY=`grep -m1 mapped bwa1/$sample.chrY_flagstat | awk '{print $1}'`
percent_mapped_chrY=`echo $tot_mapped $tot_mapped_chrY  | awk '{print $2*100/$1}'`
tot_mapped_prop_only_chrY=`grep -m1 mapped bwa1/$sample.properly_onlyChrY_flagstat | awk '{print $1}'`
percent_mapped_prop_only_chrY=`echo $tot_mapped_chrY $tot_mapped_prop_only_chrY | awk '{print $2*100/$1}'`
tot_mapped_q30_prop_only_chrY=`grep -m1 mapped bwa1/$sample.q30_properly_onlyChrY_flagstat | awk '{print $1}'`
echo -e "$sample\t$cpu\t$mem\t$tot_read\t$tot_mapped\t$tot_properly_paired\t$tot_mapped_chrY\t$percent_mapped_chrY\t$tot_mapped_prop_only_chrY\t$percent_mapped_prop_only_chrY\t$tot_mapped_q30_prop_only_chrY" >> bwa1/nb_aln.tsv
done

Pour les séquences proprement pairées et alignées uniquement sur le chromosome Y avec une mapping quality >= 30, réalignement blast pour avoir une sortie "alignement"


module load bioinfo/ncbi-blast-2.7.1+

mkdir ref_blast
ln -s ../ref_bwa1/GCA_000003205.6_Btau_5.0.1_genomic.mask.fna ref_blast/GCA_000003205.6_Btau_5.0.1_genomic.mask.fna
sbatch -J index_blast -o LOG/%x.out -e LOG/%x.err --mem=10G --wrap="makeblastdb -in ref_blast/GCA_000003205.6_Btau_5.0.1_genomic.mask.fna -dbtype nucl" 

mkdir blast

# extraction des reads
for i in bwa1/*.bam
do
echo $i
samtools view -h -f2 -q 30 $i CM001061.2 | grep -v XA:Z |samtools view -Shb - | samtools fasta - > blast/`basename $i | sed 's/\.bam/\.fasta/'`
done

# alignement
for i in blast/*.fasta
do
sbatch -J blast_`basename $i |sed 's/.fasta//'` -o LOG/%x.out -e LOG/%x.err --wrap="blastn -query $i -out `echo $i | sed 's/\.fasta/\.blast/'` -db ref_blast/GCA_000003205.6_Btau_5.0.1_genomic.mask.fna" 
done

Résultats

nb_aln.tsv

#sample_name CPU_time MEM_used nb_reads nb_read_mapped nb_read_properly_mapped nb_read_mapped_on_chrY percent_read_mapped_on_chrY nb_read_properly_mapped_only_on_chrY percent_read_properly_mapped_only_on_chrY nb_read_mapq30_properly_mapped_only_on_chrY
Emy10 00:23:56 5.62 GB 1680356 1668864 1548578 586 0.0351137 399 68.0887 16
Emy11 00:12:40 5.74 GB 635726 652803 526810 145 0.0222119 87 60 2
Emy12 00:06:45 5.29 GB 354868 362510 274570 89 0.024551 40 44.9438 0
Emy13 00:08:37 5.79 GB 450848 343499 320264 254 0.0739449 67 26.378 6
Emy14 00:22:33 6.76 GB 1401822 1081970 1043656 880 0.0813331 286 32.5 10
Emy15 00:26:29 6.48 GB 1022108 1097508 779532 509 0.0463778 264 51.8664 10
Emy16 00:21:32 6.44 GB 569334 598064 406594 158 0.0264186 70 44.3038 0
Emy17 00:10:53 6.01 GB 367904 414962 305324 164 0.0395217 86 52.439 0
Emy18 00:10:42 6.11 GB 410440 417769 375546 82 0.0196281 50 60.9756 2
Emy19 00:38:21 7.59 GB 1110116 1280700 926940 434 0.0338877 220 50.6912 6
Emy1 00:33:00 7.07 GB 1550226 1701654 1305580 385 0.022625 189 49.0909 8
Emy20 00:02:56 4.92 GB 188702 181417 174984 42 0.0231511 23 54.7619 2
Emy21 00:02:25 4.65 GB 147720 131763 126294 83 0.0629919 20 24.0964 0
Emy22 00:01:42 4.63 GB 79774 69844 63568 40 0.0572705 18 45 2
Emy23 00:12:29 6.02 GB 443806 413871 318000 220 0.0531567 105 47.7273 4
Emy24 00:25:16 6.43 GB 1042714 1040747 797926 296 0.0284411 137 46.2838 2
Emy25 00:24:36 6.73 GB 801218 783511 603822 371 0.047351 157 42.3181 8
Emy26 00:08:29 5.53 GB 390056 383760 308664 105 0.0273609 36 34.2857 0
Emy27 00:03:08 4.84 GB 201608 200321 186110 45 0.0224639 25 55.5556 0
Emy28 00:15:32 5.76 GB 1001470 1012655 890862 219 0.0216263 112 51.1416 6
Emy29 00:26:00 6.00 GB 1589610 1673287 1363234 381 0.0227696 263 69.0289 11
Emy2 00:01:00 1.38 MB 53178 50349 47714 18 0.0357505 9 50 0
Emy30 00:08:01 5.43 GB 434296 342089 302992 333 0.0973431 79 23.7237 13
Emy31 00:01:58 4.69 GB 429156 145136 141392 223 0.153649 66 29.5964 12
Emy3 00:15:46 6.27 GB 599506 635532 486484 202 0.0317844 78 38.6139 1
Emy4 00:15:49 6.21 GB 573952 650761 380622 185 0.0284283 82 44.3243 2
Emy5 00:08:51 5.70 GB 491172 511635 375054 134 0.0261905 54 40.2985 2
Emy6 00:16:26 6.05 GB 682338 804165 377532 257 0.0319586 110 42.8016 2
Emy7 00:17:35 6.12 GB 749392 835071 532060 238 0.0285006 128 53.7815 4
Emy8 00:06:12 5.37 GB 270596 246936 204308 162 0.065604 74 45.679 8
Emy9 00:28:43 6.53 GB 1483164 1436478 1328500 365 0.0254094 221 60.5479 6

report_R1.html Magnifier (1,359 Mo) Maria Bernard, 28/10/2020 16:20

report_R2.html Magnifier (1,366 Mo) Maria Bernard, 28/10/2020 16:20

nb_aln.tsv Magnifier (2,393 ko) Maria Bernard, 28/10/2020 16:36