Commit 28ee6e06 authored by Uladzislava KHAURATOVICH's avatar Uladzislava KHAURATOVICH 💬
Browse files

Update !all_species_F_M_Wsnps.md

parent 57aa065c
......@@ -5,10 +5,13 @@ We will also check for W-SNPs in parthenogenetic (asexual) species which is goin
working directory: `/nfs/scistore03/vicosgrp/ukhaurat/all_species_males_females_Wsnps`
## 1. Transform raw bam files into fasta
## 1. Transform raw bam files with RNA reads into fasta
Location of all the raw reads [Beatriz's table](https://git.ist.ac.at/bvicoso/artemia_zw_sexasex_2020/-/blob/master/UsefulData_July2021.md#rna-seq-data-head)
<details><summary>code for bam -> fastq transformation</summary>
<p>
```
bam_to_fastq
......@@ -53,10 +56,14 @@ java -jar $PICARD SamToFastq -I /archive3/group/vicosgrp/shared/artemia_RNAseq_r
java -jar $PICARD SamToFastq -I /archive3/group/vicosgrp/shared/artemia_RNAseq_reads/Atanasovsko/45191_ACTGAT_C9KRJANXX_7_20161102B_20161102.bam -F A.parAta.F.head.RNA.91_1.fq -F2 A.parAta.F.head.RNA.91_2.fq +
java -jar $PICARD SamToFastq -I /archive3/group/vicosgrp/shared/artemia_RNAseq_reads/Atanasovsko/45192_GTCCGC_C9KRJANXX_7_20161102B_20161102.bam -F A.parAta.F.head.RNA.92_1.fq -F2 A.parAta.F.head.RNA.92_2.fq +
```
</p>
</details>
Trimming
## 2. Trimming
<details><summary>code for trimming</summary>
<p>
```
module load trimmomatic
......@@ -73,11 +80,12 @@ ILLUMINACLIP:/nfs/scistore03/vicosgrp/Bioinformatics_2018/Beatriz/0-Software/Tri
SLIDINGWINDOW:4:22 MINLEN:36 LEADING:3 TRAILING:3
done
#delete tabs and enters!
/usr/bin/time -v sbatch trimming.sh
Submitted batch job 24751597
```
didn't work, needed to delete tabs and enters! fixed the loop
beforehead tried with individual files, the command worked.
```
module load trimmomatic
cd /nfs/scistore03/vicosgrp/ukhaurat/all_species_males_females_Wsnps/
......@@ -88,8 +96,17 @@ srun java -jar /nfs/scistore03/vicosgrp/Bioinformatics_2018/Beatriz/0-Software/T
/usr/bin/time -v sbatch trimming_gz.sh
Submitted batch job 24698471
```
Mapping
```
</p>
</details>
## 3. Mapping
<details><summary>code for mapping</summary>
<p>
loop never worked :( smth with naming i guess.
```
module load star
......@@ -137,10 +154,20 @@ STAR --genomeDir /nfs/scistore03/vicosgrp/ukhaurat/full_assembly_analysis/index/
STAR --genomeDir /nfs/scistore03/vicosgrp/ukhaurat/full_assembly_analysis/index/ --readFilesIn A.urmiana.M.head.RNA.77_1.trim.fq A.urmiana.M.head.RNA.77_2.trim.fq --outSAMtype BAM SortedByCoordinate --outFileNamePrefix /nfs/scistore03/vicosgrp/ukhaurat/all.species.males.females.Wsnps/A.urmiana.M.head.RNA.77_ --runThreadN 10
STAR --genomeDir /nfs/scistore03/vicosgrp/ukhaurat/full_assembly_analysis/index/ --readFilesIn A.urmiana.F.head.RNA.75_1.trim.fq A.urmiana.F.head.RNA.75_2.trim.fq --outSAMtype BAM SortedByCoordinate --outFileNamePrefix /nfs/scistore03/vicosgrp/ukhaurat/all.species.males.females.Wsnps/A.urmiana.F.head.RNA.75_ --runThreadN 10
```
</p>
</details>
__mapping statistics__
## 4. SNP calling
Genome: CHRR_integrated.fa
__I filtered for quality and coverage >10, but didn't remove multialelic sites!__
<details><summary>code for VCF file creation</summary>
<p>
### SNP calling
```
module load samtools
module load bcftools
......@@ -156,11 +183,15 @@ srun vcftools --gzvcf head.sin.kaz.urm_vcf.gz --remove-indels --maf 0.1 --max-mi
#Filter 2: remove multiallelic
#bcftools view --max-alleles 2 --exclude-types indels head_asex_raremale_cov5_filtered.vcf > head_asex_raremale_cov5_filtered2.vcf
```
sbatch snp_calling_sin.kaz.urm.sh
Submitted batch job 24779195
```
</p>
</details>
The order of individuals in the VCF file:
```
1. A.sinica.F.head.RNA.69_Aligned.sortedByCoord.out.bam
2. A.sinica.F.head.RNA.70_Aligned.sortedByCoord.out.bam
3. A.sinica.F.head.RNA.97_Aligned.sortedByCoord.out.bam
......@@ -188,3 +219,35 @@ Submitted batch job 24779195
21. A.parUrm.F.head.RNA.56_Aligned.sortedByCoord.out.bam
22. A.parUrm.F.head.RNA.57_Aligned.sortedByCoord.out.bam
```
## 5. Check for W-SNPs in unsimplified VCF
Prepare the list of W-SNP coordinates
working directory: `/nfs/scistore03/vicosgrp/ukhaurat/full_assembly_analysis/`
```
cat w_snps_perl_script_full_assembly.txt | awk '{print $1, $2}' > w_snps_coordinates.txt
mv w_snps_coordinates.txt ../all.species.males.females.Wsnps/
```
working directory: `/nfs/scistore03/vicosgrp/ukhaurat/all.species.males.females.Wsnps/`
```
grep -f w_snps_coordinates.txt head.sin.kaz.urm_filtered.vcf > w_snps_all_species.txt
#didn't work
#make a simplified vcf
#remove the filtering for NAs
cat head.sin.kaz.urm_filtered.vcf | grep -v '^##' | perl -pi -e 's/:.*?\t/\t/gi' | perl -pi -e 's/:.*//gi' | perl -pi -e 's/\t/ /gi' | perl -pi -e 's/0\/0/0/gi' | perl -pi -e 's/1\/1/2/gi' | perl -pi -e 's/0\/1/1/gi' > head.sin.kaz.urm_filtered.vcf_simple
grep -f w_snps_coordinates.txt head.sin.kaz.urm_filtered.vcf_simple > w_snps_all_species.txt
#too few
module load vcftools
vcftools --vcf head.sin.kaz.urm_filtered.vcf --snps w_snps_coordinates.txt --recode --recode-INFO-all --out w_snps_all_species.txt
#didn't work
cat w_snps_coordinates.txt | awk '{print $2}' > w_snps_positions.txt
vcftools --vcf head.sin.kaz.urm_filtered.vcf --positions w_snps_positions.txt --recode --recode-INFO-all --out w_snps_all_species.txt
#nothing works!!!!
```
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment