处理单细胞的流程

#!/bin/bash
barcode=/media/whq/282A932A2A92F3D2/WHQ/tang_analysis/HECA/barcode_96_8bp.txt
trim_script=/media/whq/282A932A2A92F3D2/WHQ/tang_analysis/HECA//trim_TSO_polyA.pl
index=/media/whq/282A932A2A92F3D2/WHQ/reference_index/ensembl_GRCH38_hisat2_index/ref.fa
gtf=/media/whq/282A932A2A92F3D2/WHQ/reference_index/ensembl_GRCH38_hisat2_index/Homo_sapiens.GRCh38.99.gtf
get_combine_matrix=/media/whq/282A932A2A92F3D2/WHQ/tang_analysis/HECA/get_combine_matrix.R
#生成测试数据
##mkdir practice_5000_fastq
##ls -d SRR*|while read id ;do cat $id/$id'_1.fastq.gz'|head -20000>practice_5000_fastq/$id'_1.fastq.gz';cat $id/$id'_2.fastq.gz'|head -20000>practice_5000_fastq/$id'_2.fastq.gz';done

#将practice_5000_fastq作为之后的fastq文件练习
#新建文件加并将文件放进去

##ls SRR*'_1.fastq.gz'|while read id ;do dir=${id%_1*};mkdir $dir ;mv $id $dir ;mv $dir'_2.fastq.gz' $dir ;done

#使用UMItools提取barcode和umi文件

ls -d SRR*|while read id ;do umi_tools extract --bc-pattern=CCCCCCCCNNNNNNNN --stdin $id/$id'_2.fastq.gz' --stdout $id/$id.R1.extracted.fq.gz  --read2-stdout  --read2-in $id/$id'_1.fastq.gz' --filter-cell-barcode --whitelist=$barcode ;done

#去除TSO和polyA


ls -d SRR*|while read id ;do perl $trim_script $id/$id.R1.extracted.fq.gz $id/$id.R1.trim.fq.gz 0;seqtk trimfq $id/$id.R1.trim.fq.gz | gzip - > $id/$id.R1.clean.fq.gz ;done

#使用hisat2跑、排序并删除sam文件

ls -d SRR*|while read id ;do read=$id/$id.R1.clean.fq.gz;hisat2 -p 15 --dta -x $index -U $read -S $id/$id'.sam';samtools sort -@ 10 --output-fmt BAM -o $id/$id'.sorted.bam'  $id/$id'.sam';  sam=$id/$id*'.sam';echo $sam;rm $sam;done




#featurecounts

ls -d SRR*|while read id;do featureCounts -T 4 -a $gtf -o $id/gene_assigned -R BAM $id/$id.sorted.bam;done

#重新使用samtools排序


ls -d SRR*|while read id ;do samtools sort -@ 10 --output-fmt BAM -o $id/$id'.sorted.bam'  $id/$id'.sorted.bam.featureCounts.bam';done

#建索引


ls -d SRR*|while read id ;do samtools index $id/$id.sorted.bam;done

#使用umitools计数

ls -d SRR*|while read id ;do umi_tools count --per-gene --gene-tag=XT --per-cell --wide-format-cell-counts -I $id/$id.sorted.bam -S $id/$id.UMI_counts.tsv;done

#获得表达矩阵！
Rscript $get_combine_matrix