samtools + bcftools

最新推荐文章于 2024-01-23 15:30:17 发布

All_Will_Be_Fine噻

最新推荐文章于 2024-01-23 15:30:17 发布

阅读量246

点赞数

分类专栏： bioinfo NGS 文章标签： NGS

本文链接：https://blog.csdn.net/jiangshandaiyou/article/details/125442400

版权

bioinfo 同时被 2 个专栏收录

65 篇文章 0 订阅

订阅专栏

NGS

19 篇文章 3 订阅

订阅专栏

samtools + bcftools处理一个donor 的800多个细胞

# conda activate ngs

/public/home/djs/software/tools_NCBI/sratoolkit.2.10.9-ubuntu64/bin/prefetch-orig.2.10.9  --option-file /public/home/djs/huiyu/20220613_scRNA_VAF/SRR_Acc_List.txt -X 1000G -p -r yes -c -O /public/home/djs/huiyu/20220613_scRNA_VAF/ -t http

ls |while read id ;do echo "~/software/tools_NCBI/sratoolkit.2.10.9-ubuntu64/bin/fasterq-dump  --split-3  ./$id/${id}.sra" >> fasterq_dump.sh; done
ParaFly -c fasterq_dump.sh -CPU 10 >out.log 2>err.log


mkdir {QC_1,QC_2,clean_data}

ls |grep ".fastq" | while read id ;do (fastqc -o ./QC_1  $id &);done

cd QC_1
mkdir multiqc
multiqc ./ -o ./multiqc

cd ../
ls |grep "1.fastq" > 1.txt
ls |grep "2.fastq" > 2.txt
paste 1.txt 2.txt > trim.txt
cat trim.txt |while read id;do a=($id) && ~/software/TrimGalore-master/trim_galore -q 25 --phred33 --stringency 3 -o ./clean_data --paired ${a[0]} ${a[1]}; done

cd clean_data/
ls |grep "report" |while read id ;do rm $id ;done
ls |grep "val.*fq" | while read id ;do (fastqc -o ../QC_2  $id &);done
cd ../QC_2/
mkdir multiqc
multiqc ./ -o ./multiqc


cd ../clean_data
ls |grep "val_1.fq" > 1.txt
ls |grep "val_2.fq" > 2.txt
paste 1.txt 2.txt  > map.txt
cat map.txt |while read id;do a=($id) && echo "/public/home/djs/software/bwa-0.7.17/bwa mem /public/home/djs/software/bwa-0.7.17/index_bwa/human_38/human_ref      ${a[0]}  ${a[1]} > ${a[0]}.sam";done >> mapping.sh
nohup ParaFly -c mapping.sh -CPU 20 >out.log 2>err.log &

# sam_to_bam && sort && index

ls |grep "fq.sam" |while read id ;do echo "samtools view -u ${id} | samtools sort -  > ${id}.sort.bam  && samtools index ${id}.sort.bam" ;done > sam_to_bam.sh

nohup ParaFly -c sam_to_bam.sh -CPU 20 >out.log 2>err.log &

# samtools mpileup && bcftools calling
ls |grep "sort.bam" | grep -v "bai" |while read id ;do echo "samtools mpileup -uf /public/home/djs/software/STAR-2.7.7a/genome_raw/GRCh38.primary_assembly.genome.fa ${id} >> ${id}.bcf && bcftools view -v snps ${id}.bcf > ${id}.vcf && java -Xmx8g -jar /public/home/djs/software/snpEff/snpEff.jar -c /public/home/djs/software/snpEff/snpEff.config  GRCh38 ${id}.vcf >> ${id}.ann.vcf";done > snv_calling.sh
nohup ParaFly -c snv_calling.sh -CPU 20 >>out.log 2>>err.log &