Calculate corrected theta from resequencing data - Revision history

Dawes001 at 15:00, 15 July 2019

2019-07-15T15:00:06Z

← Older revision		Revision as of 17:00, 15 July 2019
Line 8:		Line 8:
	#SBATCH --ntasks=1		#SBATCH --ntasks=1
	#SBATCH --nodes=1		#SBATCH --nodes=1
	#SBATCH --constraint=~~normalmem~~		#SBATCH --constraint=4gpercpu
	#SBATCH --output=output_%j.txt		#SBATCH --output=output_%j.txt
	#SBATCH --error=error_output_%j.txt		#SBATCH --error=error_output_%j.txt

Dawes001 at 14:59, 15 July 2019

2019-07-15T14:59:53Z

← Older revision		Revision as of 16:59, 15 July 2019
Line 12:		Line 12:
	#SBATCH --error=error_output_%j.txt		#SBATCH --error=error_output_%j.txt
	#SBATCH --job-name=ngstheta		#SBATCH --job-name=ngstheta
	~~#SBATCH --partition=ABGC_Research~~
	module load samtools/0.1.19		module load samtools/0.1.19
	VAR=`gunzip -c /lustre/nobackup/WUR/ABGC/shared/Pig/vars_hjm_newbuild10_2/vars-flt_$1-final.txt.gz \| cut -f8 \| head -1000000 \| sort \| uniq -c \| sed 's/^ \+//' \| sed 's/ \+/\t/' \| sort -k1 -nr \| head -1 \| cut -f2`		VAR=`gunzip -c /lustre/nobackup/WUR/ABGC/shared/Pig/vars_hjm_newbuild10_2/vars-flt_$1-final.txt.gz \| cut -f8 \| head -1000000 \| sort \| uniq -c \| sed 's/^ \+//' \| sed 's/ \+/\t/' \| sort -k1 -nr \| head -1 \| cut -f2`

Megen002 at 09:43, 1 April 2014

2014-04-01T09:43:25Z

← Older revision		Revision as of 11:43, 1 April 2014
Line 33:		Line 33:
	</source>		</source>

			Average values for Theta were then extracted with the following R-srcript:
	<source lang = 'rsplus'>		<source lang = 'rsplus'>
	files=list.files(pattern="wintheta")		files=list.files(pattern="wintheta")

Megen002 at 09:42, 1 April 2014

2014-04-01T09:42:52Z

← Older revision		Revision as of 11:42, 1 April 2014
Line 38:		Line 38:
	for (file1 in files){		for (file1 in files){
	x <- read.table(file1,header=T);		x <- read.table(file1,header=T);
	mn=mean(x$THETA_HET[x$BP>20000 & x$CHR != 'chrUN_nr']);		mn=mean(x$THETA_HET[x$BP>20000 & x$CHR != 'chrUN_nr' & x$CHR != 'Ssc10_2_X']);
	print(paste(file1,mn,sep=" "));		print(paste(file1,mn,sep=" "));
	a<- rbind(a,data.frame("file"=file1,"theta_het"=mn))		a<- rbind(a,data.frame("file"=file1,"theta_het"=mn))

Megen002 at 09:41, 1 April 2014

2014-04-01T09:41:04Z

← Older revision		Revision as of 11:41, 1 April 2014
Line 1:		Line 1:
			This procedure will estimate theta (nucleotide diversity) based on re-sequencing data. The method is describe in [http://www.biomedcentral.com/1471-2164/14/148 Esteve-Codina et al.]


	<source lang ='bash'>		<source lang ='bash'>
	#!/bin/bash		#!/bin/bash
Line 24:		Line 27:
	</source>		</source>

			The script can be submitted using <code>sbatch</code> using the following code, assuming that the names of the individuals are listed in a file called <code>individuals.txt</code>.
	<source lang='bash'>		<source lang='bash'>
	INDS=`cat ~~bams~~.txt ~~\| grep -v LB01~~`		INDS=`cat individuals.txt`
	for IND in $INDS; do sbatch nucdiv_pipeline.sh $IND; done		for IND in $INDS; do sbatch nucdiv_pipeline.sh $IND; done
	</source>		</source>

Megen002 at 10:22, 30 March 2014

2014-03-30T10:22:26Z

← Older revision		Revision as of 12:22, 30 March 2014
Line 32:		Line 32:
	files=list.files(pattern="wintheta")		files=list.files(pattern="wintheta")
	a <- data.frame("file" = character(), "theta_het" = numeric())		a <- data.frame("file" = character(), "theta_het" = numeric())
	for (file1 in files){x <- read.table(file1,header=T); mn=mean(x$THETA_HET[x$BP>20000 & x$CHR != 'chrUN_nr']); print(paste(file1,mn,sep=" "));a<- rbind(a,data.frame("file"=file1,"theta_het"=mn))}		for (file1 in files){
			x <- read.table(file1,header=T);
			mn=mean(x$THETA_HET[x$BP>20000 & x$CHR != 'chrUN_nr']);
			print(paste(file1,mn,sep=" "));
			a<- rbind(a,data.frame("file"=file1,"theta_het"=mn))
			}
	write.table(x=a,file="theta_het_results.txt")		write.table(x=a,file="theta_het_results.txt")
	</source>		</source>

Megen002 at 10:21, 30 March 2014

2014-03-30T10:21:39Z

← Older revision		Revision as of 12:21, 30 March 2014
Line 29:		Line 29:
	</source>		</source>

	<source lang = 'R'>		<source lang = 'rsplus'>
	files=list.files(pattern="wintheta")		files=list.files(pattern="wintheta")
	a <- data.frame("file" = character(), "theta_het" = numeric())		a <- data.frame("file" = character(), "theta_het" = numeric())

Megen002: Created page with " #!/bin/bash #SBATCH --time=10000 #SBATCH --mem=4000 #SBATCH --ntasks=1 #SBATCH --nodes=1 #SBATCH --constraint=normalmem #SBATCH --output=output_%j.txt #S..."

2014-03-30T10:20:47Z

Created page with "<source lang ='bash'> #!/bin/bash #SBATCH --time=10000 #SBATCH --mem=4000 #SBATCH --ntasks=1 #SBATCH --nodes=1 #SBATCH --constraint=normalmem #SBATCH --output=output_%j.txt #S..."

New page

<source lang ='bash'>
#!/bin/bash
#SBATCH --time=10000
#SBATCH --mem=4000
#SBATCH --ntasks=1
#SBATCH --nodes=1
#SBATCH --constraint=normalmem
#SBATCH --output=output_%j.txt
#SBATCH --error=error_output_%j.txt
#SBATCH --job-name=ngstheta
#SBATCH --partition=ABGC_Research
module load samtools/0.1.19
VAR=`gunzip -c /lustre/nobackup/WUR/ABGC/shared/Pig/vars_hjm_newbuild10_2/vars-flt_$1-final.txt.gz | cut -f8 | head -1000000 | sort | uniq -c | sed 's/^ \+//' | sed 's/ \+/\t/' | sort -k1 -nr | head -1 | cut -f2`
let MAX=2*VAR
echo "$1 max_depth is $MAX"
MIN=$(( $VAR / 3 ))
if [ $MIN -lt 5 ]; then MIN=4; fi

echo "$1 min_depth is $MIN"
samtools mpileup -uf /lustre/nobackup/WUR/ABGC/shared/Pig/Sscrofa_build10_2/FASTA/Ssc10_2_chromosomes.fa /lustre/nobackup/WUR/ABGC/shared/Pig/BAM_files_hjm_newbuild10_2/$1_rh.bam | bcftools view -bvcg - > $1.mig.bcf
bcftools view $1.mig.bcf | vcfutils.pl varFilter -d$MIN -D$MAX > $1.mig.vcf
awk '$6 >= 20' $1.mig.vcf > $1.miguel.vcf
samtools mpileup -Bq 20 -d 50000 /lustre/nobackup/WUR/ABGC/shared/Pig/BAM_files_hjm_newbuild10_2/$1_rh.bam | perl covXwin-v3.1.pl -v $1.miguel.vcf -w 50000 -d $MIN -m $MAX -b /lustre/nobackup/WUR/ABGC/shared/Pig/BAM_files_hjm_newbuild10_2/$1_rh.bam | ./ngs_theta -d $MIN -m $MAX > $1.wintheta
</source>

<source lang='bash'>
INDS=`cat bams.txt | grep -v LB01`
for IND in $INDS; do sbatch nucdiv_pipeline.sh $IND; done
</source>

<source lang = 'R'>
files=list.files(pattern="wintheta")
a <- data.frame("file" = character(), "theta_het" = numeric())
for (file1 in files){x <- read.table(file1,header=T); mn=mean(x$THETA_HET[x$BP>20000 & x$CHR != 'chrUN_nr']); print(paste(file1,mn,sep=" "));a<- rbind(a,data.frame("file"=file1,"theta_het"=mn))}
write.table(x=a,file="theta_het_results.txt")
</source>