\u5b98\u7f51\uff1ahttps:\/\/github.com\/gxiaolab\/scAllele<\/a><\/p>\n \u53c2\u8003\uff1ahttps:\/\/github.com\/pezmaster31\/bamtools\/issues\/135<\/a><\/p>\n split_script.py:<\/p>\n \u4e5f\u53ef\u5c1d\u8bd5\u76f4\u63a5\u7528bamtools\uff1a<\/p>\n https:\/\/github.com\/10XGenomics\/subset-bam<\/a><\/p>\n https:\/\/github.com\/10XGenomics\/subset-bam\/issues\/17<\/a> \u6b64\u5904\u6ce8\u610f -o \u7684 prefix \u4e0d\u80fd\u4e3abam\u6587\u4ef6\u7684\u540d\u5b57\uff0c\u5426\u5219\u4f1a\u5220\u9664bam\u6587\u4ef6\uff08bug?\uff09<\/p>\n 1. \u4e0b\u8f7d\u5e76\u5b89\u88c5scAllele \u5b98\u7f51\uff1ahttps:\/\/github.com\/gxiaolab\/scAllel…<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1],"tags":[],"class_list":["post-566","post","type-post","status-publish","format-standard","hentry","category-uncategorized"],"_links":{"self":[{"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/posts\/566","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/comments?post=566"}],"version-history":[{"count":18,"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/posts\/566\/revisions"}],"predecessor-version":[{"id":605,"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/posts\/566\/revisions\/605"}],"wp:attachment":[{"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/media?parent=566"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/categories?post=566"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/tags?post=566"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}conda create -n scallele\nmamba install -c bioconda -c conda-forge python=3.10 pysam samtools bamtools\n\npip install scallele<\/code><\/pre>\n2. \u6309\u7167Barcode\u5206\u5272\u6bd4\u5bf9\u540e\u5355\u7ec6\u80deBam\u6587\u4ef6(Cellranger out)<\/h3>\n
2.1 \u65b9\u6cd5\u4e00<\/h4>\n
ulimit -n 65535\nsamtools sort -t CB \/groups\/phyllodes\/home\/share\/Results\/scRNA\/cellranger_count\/FETB01-MPT-01\/outs\/possorted_genome_bam.bam > .\/CBsorted_tags.bam<\/code><\/pre>\n##### Code has not been tested on unsorted bam files, sort on barcode (CB):\n##### samtools sort -t CB unsorted.bam > sorted_tags.bam\n###\n##### INPUT: .bam file to be sorted and output directory to place split BC\n##### OUTPUT: .bam file for each unique barcode, best to make a new directory\n\n### Python 3.6.8\nimport pysam\n\n### Input varibles to set\n# file to split on\nunsplit_file = "\/path\/to\/dir\/of\/data\/sorted_tags.bam"\n# where to place output files\nout_dir = "\/path\/to\/dir\/of\/out_data\/"\n\n# variable to hold barcode index\nCB_hold = 'unset'\nitr = 0\n# read in upsplit file and loop reads by line\nsamfile = pysam.AlignmentFile( unsplit_file, "rb")\nfor read in samfile.fetch( until_eof=True):\n # barcode itr for current read\n CB_itr = read.get_tag( 'CB')\n # if change in barcode or first line; open new file \n if( CB_itr!=CB_hold or itr==0):\n # close previous split file, only if not first read in file\n if( itr!=0):\n split_file.close()\n CB_hold = CB_itr\n itr+=1\n split_file = pysam.AlignmentFile( out_dir + "CB_{}.bam".format( itr), "wb", template=samfile)\n\n # write read with same barcode to file\n split_file.write( read)\nsplit_file.close()\nsamfile.close()<\/code><\/pre>\npython split_script.py<\/code><\/pre>\n2.2 \u65b9\u6cd5\u4e8c<\/h4>\n
cp \/groups\/phyllodes\/home\/share\/Results\/scRNA\/starsolo\/FETB01-MPT-01Aligned.sortedByCoord.out.bam .\/\nbamtools split -tag CB -in FETB01-MPT-01Aligned.sortedByCoord.out.bam<\/code><\/pre>\n2.3 \u65b9\u6cd5\u4e09<\/h4>\n
conda activate subset-bam\n\ngunzip -c \/groups\/phyllodes\/home\/share\/Results\/scRNA\/cellranger_count\/FETB01-MPT-01\/outs\/filtered_feature_bc_matrix\/barcodes.tsv.gz > \/groups\/phyllodes\/home\/share\/Results\/scRNA\/cellranger_count\/FETB01-MPT-01\/outs\/filtered_feature_bc_matrix\/barcodes.tsv\nsubset-bam -b \/groups\/phyllodes\/home\/share\/Results\/scRNA\/cellranger_count\/FETB01-MPT-01\/outs\/possorted_genome_bam.bam -c \/groups\/phyllodes\/home\/share\/Results\/scRNA\/cellranger_count\/FETB01-MPT-01\/outs\/filtered_feature_bc_matrix\/barcodes.tsv.gz -o \/groups\/phyllodes\/home\/share\/Results\/scRNA\/scallele\/FETB01-MPT-01_bam<\/code><\/pre>\n2.4 \u65b9\u6cd5\u56db (\u8fd9\u4e2a\u624d\u80fd\u7528)<\/h4>\n
\nhttps:\/\/github.com\/aertslab\/single_cell_toolkit\/blob\/master\/subset_bam_per_cb.sh<\/a>
\nsubset_bam_per_cb.sh:<\/p>\n#!\/bin\/bash\n#\n# Copyright (C) 2024 - Gert Hulselmans\n#\n# Purpose:\n# Subset BAM file per cell barcode.\n\nset -e\nset -o pipefail\n\nsubset_bam_file_per_cb_chunk () {\n local bam_file="${1}";\n local barcodes_file="${2}";\n local bam_output_prefix="${3%.}";\n\n if [ ${#@} -ne 3 ] ; then\n printf 'Usage:\\n';\n printf ' subset_bam_file_per_cb_chunk bam_file barcodes_file bam_output_prefix\\n\\n';\n printf 'Arguments:\\n';\n printf ' bam_file: BAM file to subset per provided cell barcode.\\n';\n printf ' barcodes_file: File with cell barcodes to subset input BAM file.\\n';\n printf ' bam_output_prefix: Prefix used for output per CB BAM files.\\n';\n return 1;\n fi\n\n if [ "${bam_file##*.}" != "bam" ] ; then\n printf 'Error: BAM file should have ".bam" extension.\\n';\n return 1;\n fi\n\n if [ ! -f "${bam_file}" ] ; then\n printf 'Error: BAM file "%s" does not exist.\\n' "${bam_file}";\n return 1;\n fi\n\n if [ ! -f "${barcodes_file}" ] ; then\n printf 'Error: Barcodes file "%s" does not exist.\\n' "${barcodes_file}";\n return 1;\n fi\n\n if [ ! -d $(dirname "${bam_output_prefix}") ] ; then\n printf 'Error: Directory "%s" does not exist.\\n' "$(dirname "${bam_output_prefix}")";\n return 1;\n fi\n\n # First extract all requested cell barcodes from BAM file\n # before writing individual per CB BAM files.\n samtools view -@ 4 -h -D CB:"${barcodes_file}" "${bam_file}" \\\n | mawk \\\n -v bam_output_prefix="${bam_output_prefix}" \\\n -F '\\t' \\\n '\n BEGIN {\n header = "";\n CB_tag_found = 0;\n }\n {\n if ( $0 ~ \/^@\/ ) {\n # Collect SAM header.\n if ( header == "" ) {\n header = $0;\n } else {\n header = header "\\n" $0;\n }\n } else {\n CB_tag_found = 0;\n\n # Search for CB tag.\n for ( i = 12; i <= NF; i++ ) {\n if ( substr($i, 1, 5) == "CB:Z:" ) {\n CB = substr($i, 6);\n CB_tag_found = 1;\n break;\n }\n }\n\n if ( CB_tag_found == 1 ) {\n if (! (CB in CBs_found) ) {\n # Write BAM header if this is the first time this CB was seen.\n print header | "samtools view -O bam -o " bam_output_prefix "." CB ".bam -";\n CBs_found[CB] = 1;\n }\n\n # Write read to CB specific BAM file.\n print $0 | "samtools view -O bam -o " bam_output_prefix "." CB ".bam -";\n }\n }\n }\n END {\n # Close all file handles.\n for ( CB in CBs_found ) {\n close("samtools view -O bam -o " bam_output_prefix "." CB ".bam -");\n }\n }\n '\n}\n\nsubset_bam_file_per_cb () {\n local bam_file="${1}";\n local barcodes_file="${2}";\n local bam_output_prefix="${3%.}";\n\n # Number of barcodes to process in one invocation of subset_bam_file_per_cb_chunk.\n # This is used to limit the number of open file handles and the number of parallel\n # spawned samtools processes.\n local chunk_size="${4:-1000}";\n\n if [ ${#@} -lt 3 ] ; then\n printf 'Usage:\\n';\n printf ' subset_bam_file_per_cb bam_file barcodes_file bam_output_prefix [chunk_size]\\n\\n';\n printf 'Arguments:\\n';\n printf ' bam_file: BAM file to subset per provided cell barcode.\\n';\n printf ' barcodes_file: File with cell barcodes to subset input BAM file.\\n';\n printf ' bam_output_prefix: Prefix used for output per CB BAM files.\\n';\n printf ' chunk_size: Number of cell barcodes to process simultaniously.\\n'\n printf ' If more cell barcodes are provided, the input BAM file\\n';\n printf ' will be read multiple times. It is recommended to not\\n';\n printf ' set this value too high as the same number of parallel\\n';\n printf ' spawned samtools processes will be created for writing\\n';\n printf ' the per CB BAM files.\\n';\n printf ' Default: 1000.\\n';\n return 1;\n fi\n\n if [ "${bam_file##*.}" != "bam" ] ; then\n printf 'Error: BAM file should have ".bam" extension.\\n';\n return 1;\n fi\n\n if [ ! -f "${bam_file}" ] ; then\n printf 'Error: BAM file "%s" does not exist.\\n' "${bam_file}";\n return 1;\n fi\n\n if [ ! -f "${barcodes_file}" ] ; then\n printf 'Error: Barcodes file "%s" does not exist.\\n' "${barcodes_file}";\n return 1;\n fi\n\n if [ ! -d $(dirname "${bam_output_prefix}") ] ; then\n printf 'Error: Directory "%s" does not exist.\\n' "$(dirname "${bam_output_prefix}")";\n return 1;\n fi\n\n if [ $(cat ${barcodes_file} | wc -l) -le ${chunk_size} ] ; then\n # If barcodes file contains less than ${chunk_size} lines, process all barcodes at once.\n subset_bam_file_per_cb_chunk "${bam_file}" "${barcodes_file}" "${bam_output_prefix}";\n return $?;\n else\n # Split barcodes file in chunks of ${chunk_size} lines.\n local subset_bam_file_per_cb_tmp_dir=$(mktemp -t -d subset_bam_file_per_cb_XXXXXX);\n split -l ${chunk_size} -d "${barcodes_file}" ${subset_bam_file_per_cb_tmp_dir}\/barcodes_part\n\n # Subset BAM file in per chunk of barcodes.\n for barcodes_part_file in ${subset_bam_file_per_cb_tmp_dir}\/barcodes_part* ; do\n subset_bam_file_per_cb_chunk "${bam_file}" "${barcodes_part_file}" "${bam_output_prefix}";\n done\n\n # Remove temporary files.\n rm ${subset_bam_file_per_cb_tmp_dir}\/barcodes_part*;\n rmdir "${subset_bam_file_per_cb_tmp_dir}";\n fi\n}\n\nsubset_bam_file_per_cb "${@}";<\/code><\/pre>\nconda create -n subset-bam\nmamba install -c bioconda -c conda-forge mawk samtools\n\n.\/subset_bam_per_cb.sh \/groups\/phyllodes\/home\/share\/Results\/scRNA\/cellranger_count\/FETB01-MPT-01\/outs\/possorted_genome_bam.bam \/groups\/phyllodes\/home\/share\/Results\/scRNA\/cellranger_count\/FETB01-MPT-01\/outs\/filtered_feature_bc_matrix\/barcodes.tsv \/groups\/phyllodes\/home\/share\/Results\/scRNA\/scallele\/FETB01-MPT-01\n\n# \u6784\u5efa\u6279\u91cf\u8fd0\u884c\u811a\u672c\uff1a\nfor d in \/groups\/phyllodes\/home\/share\/Results\/scRNA\/cellranger_count\/*\/; do s=$(basename "$d"); echo "mkdir $s && .\/subset_bam_per_cb.sh $d\/outs\/possorted_genome_bam.bam $d\/outs\/filtered_feature_bc_matrix\/barcodes.tsv \/groups\/phyllodes\/home\/share\/Results\/scRNA\/scallele\/$s\/$s"; done > SubSetBam.sh\n<\/code><\/pre>\n3. \u5bf9\u6bcf\u4e2abam\u6587\u4ef6\u8fd0\u884cscAllele\uff1a<\/h3>\n
for f in *.bam; do base=${f%.bam}; echo "samtools index $f && scAllele -n 4 -b $f -g \/groups\/phyllodes\/home\/share\/Reference\/GRCh38_Ensembl_Ensembl104\/fasta\/genome.fa -o $base-scallele && echo $base ok"; done > run_scAllele.sh<\/code><\/pre>\n4. \u5bf9\u6bcf\u4e2a\u6837\u672c\u6587\u4ef6\u5939\u4e0b\u7684vcf\u6587\u4ef6\u7b5b\u9009PASS\u7a81\u53d8\u5e76\u6dfb\u52a0header<\/h3>\n
#!\/bin\/bash\n\n# Headers to add (with placeholder for SAMPLE)\nHEADER="##fileformat=VCFv4.2\n##FILTER=<ID=PASS,Description=\\"Variant passes all filters\\">\n##FILTER=<ID=LowQual,Description=\\"Variant does not pass one or more filtering criteria\\">\n##INFO=<ID=varType,Number=1,Type=String,Description=\\"Variant type\\">\n##INFO=<ID=varGroup,Number=1,Type=String,Description=\\"Variant Group\\">\n##INFO=<ID=varLength,Number=1,Type=Integer,Description=\\"Variant length\\">\n##INFO=<ID=varEnd,Number=1,Type=Integer,Description=\\"End position of the variant\\">\n##INFO=<ID=TandemRep,Number=1,Type=Integer,Description=\\"Count of repeat sequence adjacent to the variant\\">\n##FORMAT=<ID=DP,Number=1,Type=Integer,Description=\\"Number of reads overlapping the variant position\\">\n##FORMAT=<ID=AC,Number=R,Type=Integer,Description=\\"Number of reads containing the variant allele\\">\n##FORMAT=<ID=RC,Number=R,Type=Integer,Description=\\"Number of reads containing the reference allele\\">\n##FORMAT=<ID=AB,Number=1,Type=Float,Description=\\"Allelic balance of the variant\\">\n##FORMAT=<ID=GT,Number=1,Type=String,Description=\\"Consensus genotype\\">\n##FORMAT=<ID=GQ,Number=1,Type=Integer,Description=\\"Genotype quality\\">\n##FORMAT=<ID=RCL,Number=1,Type=Integer,Description=\\"Read Cluster id\\">\n#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT SAMPLE"\n\n# Loop through all VCF files\nfor file in *scallele.vcf; do\n echo "Processing: $file"\n\n # 1. \u63d0\u53d6cell barcode\n cell_barcode=$(echo "$file" | sed -r 's\/^.*\\.([A-Z0-9\\-]+\\-[0-9]+)\\-scallele\\.vcf$\/\\1\/')\n echo "Cell Barcode: $cell_barcode"\n\n # 2. \u521b\u5efa\u65b0\u6587\u4ef6\u540d\uff08\u6dfb\u52a0 .pass \u540e\u7f00\uff09\n new_file="${file%.vcf}.pass.vcf"\n\n # 3. \u66ff\u6362header\u4e2d\u7684SAMPLE\n header_with_cell_barcode=$(echo "$HEADER" | sed "s\/SAMPLE\/$cell_barcode\/")\n\n # 4. \u521b\u5efa\u65b0\u6587\u4ef6\uff08\u5305\u542b\u65b0header\u548c\u7b5b\u9009\u7684PASS\u884c\uff09\n # \u5148\u5199\u5165\u65b0header\n echo "$header_with_cell_barcode" > "$new_file"\n\n # \u8ffd\u52a0\u7b5b\u9009\u7684PASS\u884c\uff08\u4fdd\u7559\u539fheader\u4f1a\u5e72\u6270\uff0c\u6240\u4ee5\u8df3\u8fc7\u6240\u6709\u6ce8\u91ca\u884c\uff09\n awk '!\/^#\/ && $7 == "PASS"' "$file" >> "$new_file"\n\n echo "Created new file: $new_file (contains only PASS variants)"\ndone<\/code><\/pre>\n5. \u5bf9\u8fc7\u6ee4\u540e\u7684vcf\u8fdb\u884cAnnovar+vcf2maf<\/h3>\n
# annovar\nls *.pass.vcf | perl -ne 'chomp; my $name = $1 if ($_ =~ \/([^\\\/]+)\\.pass\\.vcf\/); print "~\/software\/annovar\/table_annovar.pl $name.pass.vcf ~\/database\/Annovar\/humandb_hg38 -buildver hg38 -out $name.pass -remove -protocol refGene,cytoBand,avsnp151,exac03,gnomad41_exome,EAS.sites.2015_08,cosmic102_unique -operation g,r,f,f,f,f,f -nastring . -vcfinput \\n"'>annovar.sh\n\n# vcf2maf<\/code><\/pre>\n6. \u5bf9\u540c\u4e00\u4e2a\u6837\u672c\u7684\u4e0d\u540c\u7ec6\u80de\u7684maf\u6587\u4ef6\u8fdb\u884c\u5408\u5e76<\/h3>\n
<\/code><\/pre>\n","protected":false},"excerpt":{"rendered":"