Demo entry 5862027

titi

   

Submitted by anonymous on Jul 25, 2016 at 18:26
Language: Perl. Code size: 20.7 kB.

## Ce programme permet l'alignement de reads issus du séquençage (454) d'ARNm de TP53. 
## L'alignement est effectué sur le génome humain (GRCh37) et les variants alpha, beta et gamma de TP53 (dont la séquence fasta est "reversée").



use warnings;
use strict;


##########################################################################################################################################################################################################################
##########################################################################################################################################################################################################################
#########################################################################################				#####################################################################################################
#########################################################################################  	     Directions		#####################################################################################################
#########################################################################################				#####################################################################################################
##########################################################################################################################################################################################################################
##########################################################################################################################################################################################################################


my $dir_genome 		= $ARGV[0]	; 	# Dossier où se trouve l'ensemble des génomes (construits via gmap build).
my $dir_analyses 	= $ARGV[1]	;	# Dossier où se trouve l'ensemble des fichiers Fastq.
my $dir_software 	= $ARGV[2]	;	# Dossier où se trouve VarScan.
my $dir_vcf 		= $ARGV[3]	;	# Dossier où les fichiers vcf sont copiés (s'ils ne sont pas vides).
my $dir_mpileup 	= $ARGV[4]	;	# Dossier où les fichiers mpileup sont copiés.
my $dir_bam	 	= $ARGV[5]	;	# Dossier où les fichiers bam sont copiés.


##########################################################################################################################################################################################################################
##########################################################################################################################################################################################################################
#########################################################################################				#####################################################################################################
#########################################################################################  	     Données		#####################################################################################################
#########################################################################################				#####################################################################################################
##########################################################################################################################################################################################################################
##########################################################################################################################################################################################################################


# Tableau contenant les noms des génomes de références préalablement construits sur lesquels vont être alignés les reads.
# Attention le premier génome n'est pas traité comme les autres : ici c'est le génome humain entier.
# Le second non plus : il correspond au variant qui sera choisi par défault si l'alignement du read est le même pour tous les variants

my @genomes = ( "GRCh37" , "p53var1rev" , "p53var3rev" , "p53var4rev" ) ;


# Table de hash dans laquelle sont enregistrées les coordonnées du premier et dernier nucléotide séquencé pour chaque amorce et ce pour chaque génome de référence.
# Ces coordonnées sont utilisées pour construire le mpileup.

my %coor = 	( 	

			"GRCh37" => 	{
			      
						"N"	=>	"17:7,578,208-7,579,526" ,
						"C"	=>	"17:7,573,983-7,578,258" 
			      
					},
					
			"p53var1rev" =>	{
			      
						"N"	=>	"gi\\|371502114\\|ref\\|NM_000546.5\\|:1,749-2,230" ,
						"C"	=>	"gi\\|371502114\\|ref\\|NM_000546.5\\|:1,346-1,799" 
			      
					},
			
			"p53var3rev" =>	{
			      
						"N"	=>	"gi\\|371502117\\|ref\\|NM_001126114.2\\|:1,882-2,363" ,
						"C"	=>	"gi\\|371502117\\|ref\\|NM_001126114.2\\|:1,346-1,932" 
			      
					},
					
			"p53var4rev" =>	{
			      
						"N"	=>	"gi\\|371502116\\|ref\\|NM_001126113.2\\|:1,809-2,290" ,
						"C"	=>	"gi\\|371502116\\|ref\\|NM_001126113.2\\|:1,346-1,859" 
			      
					}
				
		) ;


##########################################################################################################################################################################################################################
##########################################################################################################################################################################################################################
#########################################################################################				#####################################################################################################
#########################################################################################	    Déclaration		#####################################################################################################
#########################################################################################				#####################################################################################################
##########################################################################################################################################################################################################################
##########################################################################################################################################################################################################################


my ( 

$GRCh37bam ,		# Permet d'effectuer la fusion de tous les fichiers bam pour le génome GRCh37
$Mmax , 		# Permet d'enregistrer le nombre maximum de match - le nombre d'ins/del maximum pour chaque read en fonction du variant.
$choo , 		# Permet d'enregistrer sur quel variant chaque read est le mieux aligné. 
$nb , 			# Permet le calcul de $nbM et $nbO.
$nbM ,			# Nombre de nucléotides (calculé via le CIGAR) pour chaque read qui est aligné sur chaque génome de référence (sauf hg19) et ce pour chaque biopsie de chaque échantillon.
$nbO , 			# Nombre de nucléotides (CIGAR) insérés ou délétés pour chaque read pour chaque génome de référence (sauf hg19) et ce pour chaque biopsie de chaque échantillon.
@samples ,		# Tableau contenant l'ensemble des dossier-échantillons.
@biopsies ,		# Tableau contenant les différents dossiers présents dans chaque dossier-échantillon (B1 et/ou B2 et/ou P).
@fastq ,		# Tableau contenant la liste des fichiers Fastq à analyser pour chaque échantillon et chaque biopsie.
@ligne ,		# Tableau contenant la ligne courante du fichier sam traité décomposée par les tabulations.
@cig , 			# Tableau contenant le CIGAR de chaque read décomposé avec un caractère par case. 
%hash,  		# Structure complexe permettant le tri des reads en fonction du variant sur lequel l'alignement est le meilleur. 
%type 			# Table de hash avec en clé chaque type de reads (N : reads obtenus via les amorces F159 et R642, C : amorces F590 et R1045) et en valeur les fichiers bam associés.

) ;


##########################################################################################################################################################################################################################
##########################################################################################################################################################################################################################
#########################################################################################				#####################################################################################################
#########################################################################################  	     Alignement		#####################################################################################################
#########################################################################################				#####################################################################################################
##########################################################################################################################################################################################################################
##########################################################################################################################################################################################################################


@samples = `ls $dir_analyses` ;

foreach my $sample (@samples)
{
		
	chomp($sample) ;
	
	@biopsies = `ls $dir_analyses/$sample` ; 
	
	foreach my $bio (@biopsies)
	{
		
		if ($bio =~ m/Bx/)
		{
		
			chomp($bio) ;
		
			@fastq = `ls $dir_analyses/$sample/$bio/*.fastq` ;	
			
# 			Permet de créer l'entête du fichier de fusion.
			open ( OUT1 , ">$dir_analyses/$sample/$bio/rg.txt") ;
			print OUT1 "\@RG\tID:EORTC10994\tLB:MWG1\tPL:LS454\tSM:$sample$bio" ; 
			close OUT1 ;
			
			$GRCh37bam = "" ;
			
			foreach my $var(@genomes)
			{
			
				$type{'N'} = "" ;
				$type{'C'} = "" ;
				
			
				foreach my $fic (@fastq)
				{
					
					chomp($fic) ;
					$fic =~ s/\.fastq// ;
					
					if ($var eq "GRCh37" )
					{
							
						print `gmap -t 8 -D $dir_genome/$var -d $var -f samse  $fic.fastq > $fic$var.sam` ;
												
					}
					else
					{

						print `gmap --min-intronlength=15000 -t 8 -D $dir_genome/$var -d $var -f samse  $fic.fastq > $fic$var.sam` ;
						
					}
					
					print `samtools view -bS -F 4 $fic$var.sam > $fic$var.bam` ;
					print `samtools sort $fic$var.bam $fic$var` ;
					print `samtools index $fic$var.bam` ;
						
					if ($fic =~ m/159|642/ )
					{
						
						$type{'N'} .= "$fic$var.bam " ;
					
					}
					elsif ($fic =~ m/590|1045/ )
					{
						
						$type{'C'} .= "$fic$var.bam " ;
						
					}
					
					if ($var eq "GRCh37" )
					{
						
						$GRCh37bam .= "$fic$var.bam " ;
						
					}
					
				}
				
				foreach my $k (keys(%type))
				{ 
# 				
					print `samtools merge -f -rh $dir_analyses/$sample/$bio/rg.txt $dir_analyses/$sample/$bio/$sample$bio$k$var.merged.bam $type{$k}` ;
						
					print `samtools sort $dir_analyses/$sample/$bio/$sample$bio$k$var.merged.bam $dir_analyses/$sample/$bio/$sample$bio$k$var.sorted.merged` ;
					print `samtools index $dir_analyses/$sample/$bio/$sample$bio$k$var.sorted.merged.bam` ;
					
					print `samtools view -h $dir_analyses/$sample/$bio/$sample$bio$k$var.sorted.merged.bam > $dir_analyses/$sample/$bio/$sample$bio$k$var.sorted.merged.sam` ;
# 					
				}
			}
			
			
			print `samtools merge -f -rh $dir_analyses/$sample/$bio/rg.txt $dir_analyses/$sample/$bio/$sample$bio$genomes[0].merged.bam $GRCh37bam ` ;
						
			print `samtools sort $dir_analyses/$sample/$bio/$sample$bio$genomes[0].merged.bam $dir_analyses/$sample/$bio/$sample$bio$genomes[0].sorted.merged` ;
			print `samtools index $dir_analyses/$sample/$bio/$sample$bio$genomes[0].sorted.merged.bam` ;
			print `cp $dir_analyses/$sample/$bio/$sample$bio$genomes[0].sorted.merged.bam $dir_bam/$sample$bio$genomes[0].bam` ;
			print `cp $dir_analyses/$sample/$bio/$sample$bio$genomes[0].sorted.merged.bam.bai $dir_bam/$sample$bio$genomes[0].bam.bai` ;

			
			
##########################################################################################################################################################################################################################
##########################################################################################################################################################################################################################
#########################################################################################			 	######################################################################################################
#########################################################################################	Tri des sequences	######################################################################################################
#########################################################################################  en fonction du variant	######################################################################################################
#########################################################################################			 	######################################################################################################
##########################################################################################################################################################################################################################
##########################################################################################################################################################################################################################

			foreach my $k (keys(%type))
			{
			
				%hash = () ;
			
				for ( my $i = 1 ; $i <= $#genomes ; $i ++ )
				{

					
					open ( IN1 , "$dir_analyses/$sample/$bio/$sample$bio$k$genomes[$i].sorted.merged.sam") ;
							
					while (<IN1>)
					{		
							
						if ($_ =~ m/^@/)
						{
							
							$hash{'entete'}{$genomes[$i]} = "$_" ;
							
						}
						elsif ($_ ne "" )
						{
								
							$nbM = 0 ;
							$nbO = 0 ;
							
							@ligne = split ( "\t" , $_ ) ;	
							@cig = split("",$ligne[5]) ;
								
							for (my $i2 = 0 ; $i2 <= $#cig ; $i2 ++ )
							{
							
								$nb = "" ;
								
								if ($cig[$i2] eq "M")
								{
																
									for (my $i3 = $i2-1 ; $i3 >= 0 ; $i3 -- )
									{
										
											
										if ($cig[$i3] !~ /\d/)
										{
												
											last ;
												
										}
												
										$nb = "$cig[$i3]$nb" ;
											
									}
										
									$nbM += $nb ;
										
								}
							}
								
								
							for (my $i2 = 0 ; $i2 <= $#cig ; $i2 ++ )
							{
							
								$nb = "" ;
								
								if ($cig[$i2] eq "D" or $cig[$i2] eq "I")
								{
																
									for (my $i3 = $i2-1 ; $i3 >= 0 ; $i3 -- )
									{					
											
										if ($cig[$i3] !~ /\d/)
										{
											
											last ;
											
										}
										
										$nb = "$cig[$i3]$nb" ;
										
									}
											
									$nbO += $nb ;
										
								}
							}
								
								
								
							$hash{$genomes[$i]}{$ligne[0]} =	{
											
													'nbN' => $nbM - $nbO ,
													'map' => $_
											
												} ;
												
							$hash{'all'}{$ligne[0]} ++ ;		
								  
						}
					}		
							
					close IN1 ;
				}
				
				foreach my $k2 ( sort keys(%{$hash{'all'}}) )
				{
						
					$choo = "" ;
					$Mmax = 0 ;
						
					for ( my $i = 1 ; $i <= $#genomes ; $i ++ )
					{
						
						if ( exists ($hash{$genomes[$i]}{$k2}{'nbN'}) and $hash{$genomes[$i]}{$k2}{'nbN'} > $Mmax )
						{
								
							$Mmax = $hash{$genomes[$i]}{$k2}{'nbN'} ;
							$choo = $genomes[$i] ;	  
							
						}			
					}
							
					if ( $choo ne "" )
					{
						
						$hash{'print'}{$choo} .= $hash{$choo}{$k2}{'map'} ;	
						
					}
				}
		
				foreach my $var ( sort keys(%{$hash{'print'}}) )
				{

					print "$dir_analyses/$sample/$bio/$sample$bio$k$var.tri.sam\n" ;
				
					open ( OUT1 , ">$dir_analyses/$sample/$bio/$sample$bio$k$var.tri.sam") ;
					print OUT1 "$hash{'entete'}{$var}" ;
					print OUT1 "$hash{'print'}{$var}" ;
					close ( OUT1 ) ;

#########################################################################################################################################################################################################################
#########################################################################################################################################################################################################################
########################################################################################				######################################################################################################
########################################################################################   	    VarScan		######################################################################################################
########################################################################################				######################################################################################################
#########################################################################################################################################################################################################################
#########################################################################################################################################################################################################################
			
							
					print `samtools view -bS $dir_analyses/$sample/$bio/$sample$bio$k$var.tri.sam > $dir_analyses/$sample/$bio/$sample$bio$k$var.tri.bam` ;
					print `samtools sort $dir_analyses/$sample/$bio/$sample$bio$k$var.tri.bam $dir_analyses/$sample/$bio/$sample$bio$k$var.sorted.tri` ;
					print `samtools index $dir_analyses/$sample/$bio/$sample$bio$k$var.sorted.tri.bam` ;
						
					print `samtools mpileup -B -f $dir_genome/$var/$var.fasta -r $coor{$var}{$k} $dir_analyses/$sample/$bio/$sample$bio$k$var.sorted.tri.bam > $dir_analyses/$sample/$bio/$sample$bio$k$var.txt` ;
								
					print `java -Xmx4g -jar $dir_software/VarScan.v2.3.4.jar mpileup2cns $dir_analyses/$sample/$bio/$sample$bio$k$var.txt > $dir_analyses/$sample/$bio/$sample$bio$k$var.vcf --strand-filter 0 --min-coverage 5 --min-reads2 5 --min-avg-qual 30 --min-var-freq 0.045 --output-vcf 1 --variants 1`;
							
					open ( IN1 , "$dir_analyses/$sample/$bio/$sample$bio$k$var.vcf" ) ;
						
					while (<IN1>)
					{
							
						if ( $_ !~ m/^#/ )
						{
								
							print `cp $dir_analyses/$sample/$bio/$sample$bio$k$var.vcf $dir_vcf/$sample$bio$k$var.vcf` ;
			
						}	
					}
					close ( IN1 ) ;
				}
				
				print `cp $dir_analyses/$sample/$bio/$sample$bio$k$genomes[1].txt $dir_mpileup/$sample$bio$k$genomes[1].txt` ;
							
				print `samtools mpileup -B -f $dir_genome/$genomes[0]/$genomes[0].fasta -r $coor{$genomes[0]}{$k} $dir_analyses/$sample/$bio/$sample$bio$k$genomes[0].sorted.merged.bam > $dir_analyses/$sample/$bio/$sample$bio$k$genomes[0].txt` ;
					
				print `java -Xmx4g -jar $dir_software/VarScan.v2.3.4.jar mpileup2cns $dir_analyses/$sample/$bio/$sample$bio$k$genomes[0].txt > $dir_analyses/$sample/$bio/$sample$bio$k$genomes[0].vcf --strand-filter 0 --min-coverage 5 --min-reads2 5 --min-avg-qual 30  --min-var-freq 0.045 --output-vcf 1 --variants 1`;
					
				open ( IN1 , "$dir_analyses/$sample/$bio/$sample$bio$k$genomes[0].vcf" ) ;
					
				while (<IN1>)
				{
							
					if ( $_ !~ m/^#/ )
					{
								
						print `cp $dir_analyses/$sample/$bio/$sample$bio$k$genomes[0].vcf $dir_vcf/$sample$bio$k$genomes[0].vcf` ;
								
					}
				}
				
				close ( IN1 ) ;	
					
			}	
		}
	}
}

This snippet took 0.02 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).